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1.  INTRODUCTION 


The  results  reported  herein  are  the  continuation  of 
numerical  studies  of  meso-scale  phenomena  related  to  the  ef¬ 
fects  of  orography  on  momentum  and  energy  transfer  in  the 
atmosphere  and  the  interaction  of  solar  radiation  with  the 
Earth's  atmosphere. 

1.1  OROGRAPHIC  EFFECTS  ON  GLOBAL  CLIMATE 

The  scope  of  work  during  the  past  six  months  study  has 
emphasized  a  continuing  effort  to  develop  and  expand  numerical 
codes  which  may  be  used  to  understand  the  physical  processes 
which  influence  momentum  transfer.  Further  developments  on 
the  HAIFA  code  were  completed  which  enable  us  to  calculate  the 
time  dependent  effects  of  Coriolis  forces  and  turbulence  on 
two-dimensional  flow  over  mountain  ranges.  In  addition,  a 
linear  steady-state  code  based  on  an  analysis  by  Bretherton^ 
was  completed  and  several  problem  sets  investigating  the  ef¬ 
fects  of  wind  profile  and  atmospheric  stabilities  were  com¬ 
pleted.  Where  possible,  the  steady  state  results  were  compared 
with  the  time  dependent  results  of  HAIFA. 

Some  parameterization  work  was  completed  using  the 
steady-state  results  for  an  exponential  wind  profile  and  a 
constant  lapse  rate  in  the  atmosphere.  Further  work  along 
this  line  is  expected  to  be  completed  during  the  follow-on 
contract;  in  addition,  we  expect  to  develop  three-dimensional 
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codes  using  the  same  assumptions  present  in  our  two-dimensional 
codes.  This  will  allow  us  to  calculate  the  effects  of  two- 
dimensional  topography. 

1.2  RADIATIVE  TRANSFER 

The  atmospheric  radiation  code  was  completed  during 

the  past  six  months  and  calculational  -results  were  compared 

with  the  Mintz-Arakawa  parameterization  presently  used  in  the 

UCLA  global  circulation  model  (GCM) .  These  comparisons  are 

reported  in  Sections  4  and  5.  In  addition,  an  updated  descrip- 

r  2 1 

tion  of  the  code  (superceding  our  last  report1  1 )  and  a  glos¬ 
sary  of  the  input  parameters  required  to  operate  it  are  in¬ 
cluded.  Parameterization  efforts  using  the  existing  code  will 
be  the  main  emphasis  of  the  work  and  shall  represent  the  bulk, 
of  the  work  during  the  next  six  months. 

In  addition,  further  modifications  of  the  code  are 
expected  during  the  next  six  month  period.  In  particular,  we 
expect  to  study  the  effects  of  aerosols  and  the  surface 
boundary  conditions. 
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2,  MODIFICATIONS  TO  THE  HAIFA  CODE 


Two  modifications  were  made  to  the  basic  HAIFA  code 
during  the  past  six  months  efforts.  These  were  the  inclusion 
of  (1)  the  Coriolis  forces  and  (2)  an  heuristic  turbulence 
model  in  the  basic  code.  These  modifications  are  detailed  in 
this  section  of  the  report. 

2.1  CORIOLIS  TERMS  IN  2-D  MESO-SCALE  EQUATIONS 

One  of  the  major  tasks  undertaken  in  this  contract 
period  was  the  modification  of  the  HAIFA  code  to  enable  one 
to  study  the  effects  of  Coriolis  forces  on  the  flow  over 
mountain  ranges.  An  outline  of  the  basic  theory,  previously 
presented,  is  repeated  here  for  completeness. 

In  the  2-D  calculations  previously  reported,  lee  waves 
were  formed  over  the  mountain  in  a  time  interval  of  about  one 
hour,  during  which  dynamic  effects  from  the  Coriolis  force 
were  small.  (The  geostrophic  wind  used  for  the  unperturbed 
flow,  of  course,  is  strongly  influenced  by  the  Coriolis  terms.) 

If  the  mountain  range  were  more  extensive  or  if  the  fate  of 
the  waves  radiated  by  the  mountain  were  followed  to  larger 
distances  it  is  to  be  expected  that  larger  effects  from  the 
Coriolis  terms  would  be  realized. 

The  characteristic  distance  scale  for  the  Coriolis 
force  is  =  u/f,  where  u  is  a  typical  wind  speed  and  f 
is  the  Coriolis  parameter  (if  u  ^  10  m/s,  we  obtain  Lf  'v  100  km). 
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When  the  mountain  range  is  comparable  to  Lf,  an  appreciable 
modification  of  the  gravity  waves  will  result  to  form  a  com¬ 
plex  system  of  gravity-inertia  waves. 

Since  the  Coriolis  force  induces  a  turning  of  the  wind 
it  is  necessary  to  take  into  account  several  new  factors  in 
the  calculations : 

(1)  the  component  of  the  wind  parallel  to 
the  mountain; 

(2)  pressure  gradients  in  directions  par¬ 
allel  and  perpendicular  to  the  mountain 
must  be  included  to  establish  geo- 
strophic  balance  in  the  unperturbed 
flow ;  and 

(3)  the  vertical  atmosphere  structure  is 
slightly  modified  to  account  for  the 
Coriolis  contribution  to  the  hydro¬ 
static  balance  condition. 

In  the  following  formulation  we  attempt  to  parallel  the  nu¬ 
merical  treatment  of  the  HAIFA  code  as  closely  as  feasible  in 
order  to  be  able  to  compare  the  effects  of  the  Coriolis  terms 
with  those  pertaining  to  a  non-rotating  Earth. 

2.1.1  Formulation 

The  differential  equations  of  the  dry  atmosphere  are 

formulated  in  a  system  of  reference  fixed  to  a  rotating  Earth. 

r  4 1 

As  discussed  by  Thompson  we  incorporate  the  centrifugal 
terms  into  the  definition  of  the  local  gravity  to  obtain 

+  2$  x  u  =  -  ^  VP  -  kg  ,  (2.1) 
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where  U  is  the  velocity  relative  to  the  Earth's  surface, 

^  is  the  rotational  velocity  of  the  Earth,  P  is  the  pre- 
sure,  p  is  the  density  of  the  atmosphere,  g  is  the  local 
acceleration  of  gravity  which  is  assumed  to  act  in  the  verti- 

/s 

cal  direction  k  .  The  time  derivative  is  that  evaluated  fol¬ 
lowing  the  fluid  motion.  The  equations  of  mass  and  energy 
conservation  are  not  affected  by  the  Coriolis  force. 

Neglecting  the  curvature  of  the  Earth's  surface  (and 
the  resulting  centrifugal  terms  associated  with  the  relative 
velocity),  Eq.  (2.1)  can  be  resolved  into  components.  We 
choose  a  Cartesian  coordinate  system  in  which  the  x-axis  lies 
in  the  surface  and  forms  an  angle  <J>  with  the  eastward  direc¬ 
tion,  the  y-axis  lies  in  the  surface  at  the  same  angle  <}> 
with  the  northward  direction,  and  the  z-axis  is  perpendicular, 
positive  upward.  We  shall  subsequently  assume  that  the  x-axis 
is  perpendicular  to  the  2-D  mountain  range  which  is  oriented 
at  the  angle  <f>  with  the  northward  direction  (see  figure 
below) .  Denoting  x,  y,  and  z  components  of  the  velocity  by 
u,  v,  w,  the  component  equations  are: 
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5T  -  2ft  Sine  v  +  2ft  cose  cosd)  w  =  -  -  —  , 
at  r  p  9x  ' 


dv 

dt 


+ 


2ft  sine  u  -  2ft  cose 


sin<j)  w  = 


I  9P 
P  9y 


9 


dw 

dt 


2ft  cose 


(u  cosij)  -  v  sin(f>) 


I  9P 
p  3z 


g 


(2.2) 


where  ft  is  the  magnitude  of  the  rotational  velocity,  6  is 
the  latitude  of  the  position,  and  the  time  derivatives  are 
those  formed  following  the  fluid  motion. 

In  the  absence  cf  the  mountain  varrier  we  consider 
the  atmosphere  to  be  in  geostrophic  balance;  the  motion  is 
unaccelerated,  the  vertical  velocity  component  w  is  zero, 
and  the  pressure  gradients  are  just  balanced  by  the  Coriolis 
and  gravity  terms.  Denoting  the  geostrophic  state  by  sub¬ 
scripts  g  we  obtain 


•2ft  sine  v  = 
g 


—  ^  , 

Pg  ax  ' 


2ft  sinS  u  = 

g 


i  9P 

l _ 2 

pg  ' 


•2ft  cose  (u  cos*  -  V  sind>)  +  g  = 

g  g 


Pg  9Z 


(2.3) 


Since  the  velocity  of  the  unperturbed  state  is  taken  to  be 
independent  of  x  and  y  ,  Eq.  (2.3)  shows  that  the  pres¬ 
sure  at  most  need  depend  linearly  on  x  and  y 

When  the  mountain  is  present  the  pressure,  density 
and  velocity  are  perturbed  from  their  geostrophic  values.  It 
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will  be  convenient  to  introduce  the  deviation  P'  o£  the 
pressure  from  the  geostrophic  value: 


Wc  also  introduce  the  Boussinesq  approximation,  in  which  the 
departure  of  density  from  the  geostrophic  value  is  taken  into 
account  only  in  the  gravity  term  of  Eq.  (2.3).  In  all  other 
terms,  we  use  the  geostrophic  density: 


du 


cTt  “  sine  (v-v  )  +  2ft  cose  cos<}>  w  =  -  i— 


9P ' 
9x 


dv 

cTF  +  2^  sine  (u— u  )  -  2 ft  cose  sind)  w  =  -  i. 

P , 


9P' 

9y 


dw 

dt  “  cos6  [  (u-ug)cos<Ji  -  (v-v  )  sin<|>]  =  -  2Z1  +  A  _  P_  \ 

9  pg  8z  \  '  Pg  /9  • 

(2.5) 

These  equations  constitute  the  Boussinesq  approximation  for 
e  equations  of  motion  when  the  geostrophic  flow,  assumed 
steady  and  independent  of  x  and  y  ,  is  perturbed,  clearly, 
they  are  only  approximately  satisfied  in  a  local  region,  since 
Ug  and  vg  are  not  constant,  on  the  synoptic  scale. 

We  now  consider  the  special  case  in  which  the  initial 
and  boundary  equations  are  independent  of  the  y-coordinate 
corresponding  to  a  uniform  but  obliquely  incident  wind  en¬ 
countering  a  two-dimensional  ridge,  the  topography  of  which 
is  independent  of  y  .  m  this  case,  the  initial  and  bound¬ 
ary  conditions  and  the  equations  depend  only  on  the  coordi- 

r.ates  x  and  z  .  The  resulting  2-D  equations  derived  from 
Eq.  (2.5)  are: 
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3t  +  u3jc  +  w3"z  ”  sin0  (v-v  )  +  2ft  cos0  cos*  w 


1  3p ' 


3v  .  3v  ,  3v  __ 

Jt  ulTx  +  WTS'  +  2Q  sine  (u-u  )  -  2ft  cose  sin*  w  =  0  , 


3w  9w  3w 

at  +  u3^  +  wa7  +  cose  [  (u-ug)  cos*  -  (v-vg)sin<f>]  = 


1_  _3PJ 
pg  3z 


(2.6) 


Equations  (2.6)  are  to  be  supplemented  with  the  equa¬ 
tions  of  incompressibility  and  the  temperature  equation „  The 
equations  are  unchanged  from  those  in  Ref.  2;  they  appear  in 
that  document  as  Eqs.  (2.3)  and  (2.10). 

There  are  two  major  modifications  of  the  equations 
which  have  resulted  from  the  treatment  of  the  rotation  of  the 
Earth: 

(1)  The  component  of  the  wind  parallel  to  the  range 
influences  both  the  x  and  z  momentum  equations  through 
the  Coriolis  terms.  The  y  momentum  component  equation  is 
only  weakly  coupled  to  the  others  and  may  be  solved  in  a 
similar  way  to  the  temperature  equation. 

(2)  Additional  terms  from  the  Coriolis  force  enter 
the  momentum  equations  giving  rise  to  new  terms  in  the  vorti- 
city  equation. 

2.1.2  Difference  Equations 

The  difference  equation  formulation  corresponding  to 
Eq.  (2.6)  and  supplementary  equations  can  be  chosen  to  parallel 
that  of  the  HAIFA  code.  The  equation  for  the  y-component  of 
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vorticity  is  obtained  by  cross-differentiation  of  the  x  and 
z  components  of  the  momentum,  thereby  eliminating  the  pres¬ 
sure  from  the  equations  entirely.  The  resulting  equations  are: 

If  +  ul£  *  v4£  -  2f!  Eine(lf  -  3T2-)  -  2!1  cose  sin*!!  =  p"  If  - 
3v  '  3v  8v 

3t  +  u3x  +  wa7  +  sine  *u“ug)  -  2n  cose  sin<j>  w  =  0  ,  (2.7 

where  the  y-component  of  vorticity,  n,  is  defined  by 


3u 
Tz  " 


A  quantity,  v',  is  now  introduced  which  is  the  difference 
between  the  y-component  of  wind  velocity  and  the  y-component 
of  the  geostrophic  wind, 

V  =  v  -  vg  . 

The  system  of  equations  that  HAIFA  solves  now  becomes 


v2ijj 

=  n  :  u  -  H 

co  ko 

II 

m 

i 

8T ' 

,  a(uT')  .  a(wT') 

3T 

o  .  _ 

at 

ax  az 

“  WaT"  +  F 

• 

r 

in  + 

at 

e(un)  ,  a(wn) 
ax  az 

-  20  cosS  sind)  | —  =  2—  i£ 
y  3x  p  8x 
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3v'  8 (UV1 )  3  (Wv 1 ) 

3t  3x  3z 


,3.v 

w^— 2. 
3z 


(2.8) 


-  2ft  sin0  (u-u  )  +  2ft  cos0  sin<j)  w  =  0 
y 

The  difference  approximations  for  the  above  equations  follows 
that  of  the  basic  HAIFA  formulation.  The  time-dependent 
equations  are  solved  in  an  explicit  two-level  formulation. 

The  advection  terms  are  obtained  by  using  a  high  order  con¬ 
servation  scheme  in  which  the  two  directions  are  integrated 
by  the  splitting  technique.  Additional  terms  are  all  centered 
in  space  from  quantities  available  at  the  current  time  cycle. 
The  additional  equation  for  v'  is  similar  to  the  temperature 
equation  in  structure.  The  vorticity  equation  is  modified  by 
gradients  in  the  velocity  perturbation.  Changes  in  v'  are 
due  to  both  forcing  terms  which  are  the  result  of  the  Coriolis 
force  actions  on  the  u-  and  w-components  of  velocity  and  ad¬ 
vection  of  the  Vg-component  of  geostrophic  velocity  in  the 
vertical  direction.  Depending  on  the  geostrophic  wind  struc¬ 
ture  the  advection  term  could  dominate  the  change  in  the  v' 
parameter. 

2.1.3  Initial  Conditions 

The  initial  conditions  for  HAIFA  modified  to  include 
Coriolis  forces  are  the  same  as  basic  HAIFA  except  that  the 
angle  between  the  mountain  range  and  the  north-south  direction 
and  the  angle  between  the  geostrophic  wind  and  the  mountain 
range  must  be  specified. 

The  figure  in  section  2.1.1  defines  the  parameters  used 
to  describe  the  orientation  of  the  geostrophic  wind  and  the 
mountain  range.  The  angles,  4>  and  0,  are  both  specified  in 
degrees;  <p  is  the  angle  that  the  mountain  range  makes  with  the 
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north-south  coordinate,  and  3  is  the  angle  the  geostrophic 
wind  makes  with  respect  to  the  direction  normal  to  the  mountain 
range.  The  x-  and  y--axis  in  the  code  corresponds  to  the  direc¬ 
tion  along  the  mountain  range  and  normal  to  the  mountain  range, 
respectively. 

The  other  parameter  required  in  the  Coriolis  terms  is 
the  latitude,  0,  which  is  specified  in  the  usual  manner,  in 
degrees  from  the  equator. 

The  initial  geostrophic  wind,  u  ,  is  simply  the  wind 
•  9 

which  would  be  present  without  the  mountain.  This  wind  is 

assumed  not  to  var^  in  direction  to  the  vertical.  Thus,  it  is 

similar  to  the  initial  u  velocity  in  basic  HAIFA  except  that 

the  wind  is  not  necessarily  normal  to  the  mountain.  Noting 

that 


ug(z)  =  u  ( z )  x  +  v(z)  y  (2.9) 

and  referring  to  the  figure  above  we  can  calculate  the  initial 
u-  and  v-velocities  by 


u(z)  =  cos 3  u  (z) 
g  ' 

v(z)  =  sin3  u  (z) 

The  gradient  of  the  u  -velocity  in  the  vertical  is  simply 

y 


9  u  9  u  ( z ) 

li2  =  sin6  — ! — 


(2.10) 


2.1.4  Test  Problem 

Several  preliminary  runs  were  made  using  HAIFA  with 
Coriolis  forces.  The  first  runs  were  made  with  the  angular 
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velocity,  ft,  set  to  zero  and  with  the  initial  wind  perpendicular 
to  the  mountains  :.n  an  attempt  to  duplicate  the  double  wave 
problem  run  with  the  basic  HAIFA  code. 

After  verifying  that  HAIFA  with  Coriolis  forces  dupli¬ 
cates  the  standard  double  wave  problem,  additional  runs  were 
made  with  the  angular  velocity,  ft,  set  to  its  nominal  value, 
and  with  the  wind  perpendicular  to  the  mountain  range  (i.e., 
3=0).  The  results  of  this  run  showed  virtually  no  effects 
of  the  Coriolis  terms  incorporated  in  HAIFA.  This  result  is 
not  unexpected  when  we  recall  that  the  mountain  range  in  the 

double  wave  problem  is  only  625  meters  high  and  4500  meters 
long. 

However,  problems  simulating  much  larger  mountains, 
for  example,  the  Andes  or  the  Sierras,  should  be  significantly 
affected  by  the  Coriolis  forces,  and  the  resulting  lee  waves 
produced  would  be  a  system  of  both  inertia  and  gravity  waves. 
These  complex  wave  systems  could  appreciably  change  the  drag 
from  what  one  would  expect  considering  only  gravity  waves. 

The  test  problems  were  not  continued  pending  completion 
of  a  linear  steady  state  code  described  in  Section  3  of  this 
report.  The  primary  reason  for  this  was  the  desirability  of 
determining  the  relationship  between  steady  state  and  transient 
results  of  the  vertical  flux  of  horizontal  momentum.  In 
particular,  the  cyclic  behavior  of  the  Reynolds  stress  as  a 
function  of  time  previously  reported  must  be  related  to  the 
steady  state  value  and  further  test  calculations  were  delayed 
in  order  to  attempt  to  determine  this  relationship. 
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2.2  HEURISTIC  NUMERICAL  MODEL  OF  TURBULENCE 

Xn  incompressible  flow,  whether  laminar  or  turbulent, 
the  equations  of  momentum  and  continuity,  along  with  the 
boundary  and  initial  conditions,  suffice  to  establish  com¬ 
pletely  the  exact  fluid  motion.  If  the  flow  happens  to  be 
turbulent,  however,  the  motion  involves  such  small  and  rapid 
changes  that,  although  it  is  in  principle  determinate,  its 
actual  calculation  would  impose  an  overwhelming  computational 
burden  and,  in  addition,  the  detailed  initial  conditions  are 
not  known. 

The  usual  way  around  this  difficulty  is  to  average 
the  equations  of  momentum  and  continuity  to  obtain  mean  flow 
quantities  which  are  smooth.  This  results,  of  course,  in  an 
enormous  simplification.  Unfortunately,  it  also  involves  a 
significant  and  irretrievable  loss  of  essential  information. 
Consequently,  owing  to  the  presence  of  the  unknown  Reynolds 
stresses  which  are  created  by  this  averaging  process,  the 
averaged  equations  of  momentum  and  continuity  do  not  in  them¬ 
selves  comprise  a  determinate  set.  Additional  relations  are 
required  to  fix  the  unknown  Reynolds  stresses.  Equations  can 
also  be  derived  for  the  fluctuations  and  for  averages  of  pro¬ 
ducts  of  them.  However,  these  equations  form  a  coupled  sys¬ 
tem  involving  ever  higher  order  variances.  No  first -principles 
method  for  terminating  this  system  of  equations  is  known. 

A  plausible  heuristic  approach  has  been  developed  by 
Gawain  and  Pritchett.  1  5]  Since  the  necessary  supplementary 
relations  cannot  be  established  from  the  original  equations 
by  an  analytic  procedure,  Gawain  and  Pritchett  closed  the  sys¬ 
tem  with  the  addition  of  empirical  hypotheses.  Their  philoso¬ 
phy  and  rationale  is  quoted  below: 

"From  another  viewpoint,  it  may  be  stated  that  the 
averaged  equations  of  motion  show  the  effect  of  the  Reynolds 
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stresses  upon  the  mean  flow.  However,  the  reciprocal  effect 
of  the  mean  flow  upon  the  Reynolds  stresses  is  lost  in  the 
averaging  process.  Hence  some  adequate  hypothesis  must  be 
found  for  representing  this  relation,  at  least  approximately." 

"For  this  purpose,  a  heuristic  approach  which  seems 
plausible  is  to  postulate  a  relation  between  the  Reynolds 
stresses  and  the  mean  flow  which  is  analogous  to  the  relation 
that  is  known  to  govern  the  viscous  stresses.  The  analogue 
of  the  ordinary  molecular  kinematic  viscosity  is  the  so-called 
eddy  kinematic  viscosity.'  The  problem  becomes,  therefore,  to 
determine  empirically  the  general  law  which  governs  this  mean 
effective  eddy  viscosity  at  every  space/time  point  in  the 
flow  field." 

"The  eddy  viscosity  presumably  depends  on  a  number  of 
variables,  one  of  the  most  important  of  which  is  the  local 
kinetic  energy  of  turbulence.  Therefore,  it  becomes  necessary 
to  find  the  space/time  distribution  of  the  turbulent  energy. 
Fortunately,  the  governing  energy  equation  can  be  deduced 
rigorously  from  the  original  equations  of  motion.  However, 
the  energy  equation  itself  introduces  two  additional  unknowns 
which  can  only  be  approximated  in  the  same  heuristic  and 
empirical  fashion  as  the  eddy  viscosity  itself.  The  addition¬ 
al  unknowns  are  the  rate  of  dissipation  of  turbulent  energy 
into  heat,  and  the  rate  of  turbulent  diffusion  of  energy." 

"Theory  and  experiment  both  show  that  the  eddy  viscos¬ 
ity,  and  the  dissipation  and  diffusion  functions  as  well, 
depend  not  only  on  the  turbulent  energy  itself,  but  also  on  a 
local  length  scale  parameter  which  can  be  associated  with 
each  space/time  point  in  the  flow  field.  Von  Karman  was  per¬ 
haps  the  first  to  point  out  how  a  physically  meaningful 
characteristic  length  can  be  defined  in  terms  of  local  space 
derivatives  of  the  mean  velocity  at  any  point  in  the  flow. 
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In  the  present  paper,  the  original  approach  of  von  Karman  is 
further  developed  and  refined.  It  now  takes  into  account  not 
only  the  velocity  derivatives  at  the  designated  point  itself, 
but  also  the  values  in  the  general  vicinity  of  the  point." 

"By  employing  dimensional  analysis,  and  by  applying 
the  available  experimental  data,  we  finally  obtain  three  empir¬ 
ical  expressions  which  determine  to  a  reasonable  approximation 
the  eddy  viscosity,  the  heat  dissipation  and  the  turbulent 
diffusion,  respectively.  These  expressions  also  involve  the 
turbulent  energy,  the  local  length  parameter,  and  the  distance 
to  the  nearest  fixed  wall  (if  any) .  Of  course,  these  empiri¬ 
cal  expressions  are  amenable  to  further  investigation  and 
development." 

"In  this  way  a  single  consistent  and  determinate  set 
of  equations  is  established  which  applies  in  principle  to  any 
incompressible  turbulent  flow  field.  Only  the  boundarv  condi¬ 
tions  differ  for  each  specific  application." 

2.2.1  Formulation 

The  formulation  of  equations  as  developed  by  Gawain 
and  Pritchett  must  be  modified  for  flow  of  an  incompressible 
fluid  to  flow  of  a  fluid  characterized  by  the  Boussinesq 
approximation.  Starting  with  the  momentum  and  continuity 
equations  shown  below,  the  derivation  of  equations  follows. 

+  37  =  0  (continuity)  (2.11) 


3u  ,  3u  '  3u  1  ‘  3p  ,  n  ...  „  . 
^_  +  u_  +  w__.  _  _£  +  s/  •  (Kv  Vu) 

o 


3w  ,  ?w  .  9w  1  3p  ,  „ 
=  ~  - -  TT*-  +  V  • 


(2.12) 


(momentum) 


TFt  +  +  w3¥ 


fi  +  V- (Kv  Vw)  -  f.  g 


(2.13) 
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Here  Kv  is  the  molecular  viscosity.  The  velocity  components, 
the  pressure,  and  density  can  be  separated  into  mean  and  fluc¬ 
tuating  parts, 

u  =  u  +  u' 


w  =  w  +  w' 

P  =  P  +  P' 

P  =  "p  +  P' 


These  are  to  be  inserted  into  the  equations  of  motion  and 
ensemble  averaged.  Eqs.  (2.  )  and  (2.  )  give  the  mean 

momentum  equations. 


H  +  “-i  +  k  <^>  +  "i  +  h 


and 


=  -  —  +  V *K  ,  Vu 

.  pQ  3x  v 


3w  .  — 3w  .  3  <  ■  ■  %  .  ■— 3w  .  3  ,  1 2  \ 
ITt  +  U 3x  +  3x  (u  W  >  +  W37  +  37  (w  } 


=  -  —  +  v*K  ,  Vu  +  S—  g 

pQ  3z  v  PQ 


(2.14) 


(2.15) 


Following  Gawain  and  Pritchett,  the  Reynolds  stresses 
are  postulated  to  be  related  to  the  strain  rates  of  mean  flow 
through 
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_  _  1  i 


.  .3u^  3.u . 

-u:u’  =  -  ^  u/u.'fi,  .  +  e  ir-i  +  — i 
1  vl  3  k  k  i]  3x^  3X_. 


(2.16) 


where  6i;.  -  0  for  i^j ,  and  =1  otherwise;  e  is  called  the 
eddy  kinematic  viscosity. 

The  postulate  of  Eq.  (2.16)  provides  four  necessary 
relationships , 

-u'  2  =  -  j  (u'  2  +  ^7T  +  +  2e  || 


-u 1 w'  =  e 


3w  3u 
3x  3z 


-w.'  2 


=  -  i  (u1 2  +  v' 2  +  w' 2)  +  2c 


-w'u'  =  -  u ' w ' 


(2.17) 


Applying  the  postulate  of  Eq.  (2.16)  to  Eq.  (2.14)  and 
substituting  the  kinematic  pressure 


$  =  P  “  j  u' 2  +  v' 2  +  w' 2 


one  has 


3u  +  tt3u  j.  “3u 

at  +  u53f  + 


-  (t+Kv)72u 


+  9  3e  3u  3e  /3w  ,  3u\ 

+  2  a*  a*  +  a?  (a?  +  si-]  - 


and  for  Eq.  (2.15)  one  obtains 
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9w  —aw  — 3w  1  3$ 
3t  u3x  w3z  ~  p  3z 

r  rv 


(e+K  ).V2w  + 
v 


3e  /3w  3u\  ?  3_t  3w 

ax  yax  dz)  a z  az 


(2.19) 


Equations  (2.18)  and  (2.19)  are  reformulated  in  terms 
of  the  vorticity  n  and  the  stream  function  \p,  which  satisfy 

the  relations  n  =  and  ,V2\|;  =  n-  The  vorticity  equa¬ 

tion  then  has  the  form 


Dn  = 
Dt 


(e+Kv) V2n 


2 _  IP  +  1*1$.  -  lit)  lilt  -  lit  \ 

pQ  Ex  \Sz2  ax2/  ^a.z 2  3xz  J 


..  j  a 2 e  a2ij)  a_t  an  ,  at  an_  a_e  L2  a^ 

3x3z  3x3z  3z  az  3x  3x  3x  3z 


(2.20) 


The  energy  equation  is  formed  by  multiplying  the 
momentum  equation  by  velocity.  As  before,  the  resultant  equa¬ 
tion  is  rewritten  with  mean  and  fluctuating  terms,  and  the 
averaged  energy  equation  is  subtracted  to  yield  the  turbulent 
energy  equation.  It  is  convenient  to  express  the  result  in 
Cartesian  tensor  notation. 
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<(> '  is  the  perturbation  of  kinematic  pressure,  pressure/ 
density.  The  terms  on  the  right-hand  side  of  the  energy  equa¬ 
tion  represent,  respectively,  turbulent  energy  production  cor¬ 
responding  to  the  work  done  by  the  mean  flow  against  the 
Reynolds  stresses,  dissipation  of  turbulent  energy  to  heat, 
turbulent  diffusion  of  energy,  and  molecular  diffusion.  For 
problems  at  high  Reynolds  number  the  last  term  is  vanishingly 
small;  it  will  be  ignored  hereafter. 


Denoting  by  E  the  turbulent  kinetic  energy. 


1 

2 


(u' 2  +  v' 2  +  w' 2) 


(2.22) 


and  applying  the  postulate  of  Eq.  (2.16)  to  Eq.  (2.21),  the 
turbulent  energy  equation  becomes 
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At  this  point,  neither  the  vorticity  equation  nor  the 

turbulent  energy  equation  are  closed.  It  should  also  be  noted 

here  that  the  development  of  the  turbulent  energy  equation  (2.23) 

as  described  by  Gawain'  and  Pritchett  neglects  the  temperature 

stratification  term  -  p  *  u^  -2  on  the  right-hand  side  of  the 

*o 

equation,  which  for  our  applications  can  be  a  significant  effect. 
For  the  initial  development  of  the  turbulent  scheme,  however, 
this  term  was  also  ignored  as  it  requires  that  the  fluctuating 
density  component  be  .described  in  some  heuristic  and  empirical 
fashion,  i.e.,  that  p ' u!  5  e'  ■!£•.  The  value  of  e'  remains  to 

*  1  o  Z 

be  determined  from  experimental  data  in  much  the  same  manner  as 
the  kinematic  eddy  viscosity  e  itself  needs  to  be  determined. 

As  work  on  the  turbulence  scheme  itself  progresses,  investiga¬ 
tions  of  this  term  will  be  made  through  literature  searches  in 
order  to  incorporate  its  effect  in  mountain  wave  problems. 

To  enable  the  kinematic  eddy  viscosity  e  to  be  deter¬ 
mined,  Gawain  and  Pritchett  postulated  the  formulation 


e  =  aX  /2E  .  (2.24) 

where  a  is  a  dimensionless,  slowly  varying  universal  function 
not  predictable  from  theory,  but  estimable  from  experimental 
data  and  X  is  a  length  scale  of  turbulence  in  the  vicinity  of 
a  point.  Thus,  there  is  a  X  associated  with  every  point  of 
the  flow  field.  It  was  hypothesized  that  the  definition  of  X 
in  the  vicinity  of  an  arbitrary  point  should  depend  only  on 
the  mean  flow  conditions  in  a  finite  region  surrounding  that 
point.  By  use  of  a  weighting  function  which  falls  off  rapidly 
with  increasing  separation,  dependence  on  all  points  in  the 
flow  field  can  be  avoided. 
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The  weighting  function  chosen  was 


w(x)  = 


/  '  Ax»'Ax\ 

\  A  lx  )  / 

M  M) 


all 

space 

A  strain  rate  tensor  was  defined. 


(2.25) 


r  =  +  lii 

ij  9z  3x 


(2.26) 


Note  that  =  0  for  the  incompressible  case  from  continuity. 

Next  a  generalized  strain  rate,  ft2,  an  a  generalized  strain 
rate  gradient,  ft'2,  we re  defined  as  follows: 


ft2  =  4  r.  .r. . 

2  13  13 


ft 


I  2  _  /  3ft\  7  3ft \ 
-  \3xj 


(2.27) 


(2.28) 


It  was  also  found  useful  to  define 


=  t  t)  (Hj) 


(2.29) 


A2  can  be  defined  in  terms  of  ft2  and  (ftft')2 


A2 (x)  =  I2 (x)/J2 (x) 


(2.30) 


where 
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and 


I2(x)  = 


/ 


A  -f 


w(x,x'  )n4  (x)dv' 


all 

space 


J2(x)  = 


/ 


/v  ,-K  ->■ 


w  (x,x' )  (Aft1  (x) )  2dv' 


all 

space 


expL- 


w(x,x')  = 


(x-x1 ) »  (X-X1  )  \ 
A2  (x)'  / 


A  (X) 


/  •»(-  "SM 


all 

space 


(2.31) 


(2.32) 


(2.33) 


The  last  term  in  the  turbulent  energy  equation  (Eq.  2.23) 
is  not  in  a  form  amenable  to  calculation.  It  was  postulated 
to  be  expressible  in  the  form 


u: 


/utu'.  \ 


ye 


(2.34) 


Dimensional  considerations  suggested  to  Gawain  and  Pritchett 
that  the  dissipation  of  turbulent  energy  into  heat  could  be 
expressed  in  the  form 


•  Kv/^j  8uif  V2E> 

H  2  \3xf  3x.  /  ,2 

1  3  An 


(2.35) 


which  amounts  to  a  definition  of  a  dissipation  length  AD.  To 
define  Ap,  Gawain  and  Pritchett  turned  to  the  results  of  ex¬ 
perimental  studies  at  high  Reynolds  number,  where  heat  dis¬ 
sipation  effect  tended  to  become  independent  of  Reynolds 
number.  Two  lengths  Lj  and  L2  were  defined, 
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(2.36) 

(2.37) 

and  the  relation 

JjjLj 

-p-  =  3  (2.38) 

D 

was  postulated.  The  energy  dissipation  term  becomes 

=  8  (2E)  7/6jV3  .  (2.39) 

The  complete  turbulent  energy  equation,  including  heuristic 
substitutes,  is  then  written  as  follows: 

9E  9  —  1 

9t  +  9x  (uE)  +  dz  (™E)  =  /ZE  ft2  -  8(2E)7/eJl/3 

+  4  («*  n*  H) +  h  («r*  ^  If)  (2-4o> 

To  complete  the  formulation,  it  is  necessary  to  specify 
the  computation  of  the  three  dimensionless  coefficients,  a,  8/ 
and  y»  Gawain  and  Pritchett  used  the  following  expressions, 
based  on  experimental  data: 

a  =  0.065  (1  +  exp.[-  (y/A-1)  2]  } 

1/3  =  3.7  U  +  exp  [-(y/A-1)2]}  (2.41) 

Y  =  1.4  -  0.4  exp  [- (y/A-1)  2] 

where  y  is  the  distance  to  the  nearest  fixed  boundary. 
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Thus,  the  modifications  to  the  usual  incompressible 
formulation  required  to  incorporate  the  heuristic  turbulence 
scheme  of  Gawain  and  Pritchett  are 

(1)  the  revised  vorticity  equation  [Eq.  (2.20)  ],  ani 

(2)  the  inclusion  of  an  heuristic  equation  for 
turbulent  energy  (Eq.  2.21). 

These  modifications  have  been  carried  out,  and  the 
results  are  detailed  in  the  following  sections. 

2.2.2  Numerics 

The  modifications  to  the  vorticity  ecruation  were  coded 
and  incorporated  using  typical  finite  difference  formulations, 
and  treating  the  sum  (Kv+e)  as  a  total  diffusion  term. 

The  solution  of  the  turbulent  energy  equation  was 
carried  out  in  a  new  code  subpackage,  subroutine  TURB.  The 
flow  logic  of  this  subroutine  is  displayed  in  Figure  2.1. 

The  turbulent  energy  code  was  debugged  exercised  on 
some  test  problems.  It  became  apparent  that  a  major  calcula- 
tional  burden  was  imposed  by  the  development  of  the  I2  and  J2 
terms,  and  the  associated  weighting  terms  required  for  each 
grid  point.  This  is  despite  the  fact  that  the  potentially  ex¬ 
pensive  exponential"  evaluations  can  be  reduced  to  a  one  pass 
computation,  with  tabular  evaluation  thereafter  each  cvcle. 

The  cost  in  computer  time  to  simply  form  all  the  I2  and  J2  terms, 
when  only  the  nearest  36  cells  are  used  in  the  weights,  is  about 
12  seconds.  The  cost  of  a  complete  calculational  cycle  for  the 
turbulent  formulation  approached  7  seconds  in  this  case,  nearly 
10  times  the  cost  of  the  standard  HAIFA  solution.  In  one  test, 
the  problem  characteristics  were  such  that  the  weights  did  not 
fall  off  sufficiently  rapidly,  and  it  was  required  to  include 
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SUBROUTINE  TURB 


foem  If  >  If- 


FORM  1 2  AND 

)  2 


f IS  A  NEW  > 
ESTIMATE  M 
“STIE3IRED7  > 


CALL  SECOND  ORDER 
ADVECTION  ROUTINE 


EVALUATE  EXPONENTIALS 
I  Ax2  I 

exp  -  ^  FOR 


FORM  UPDATED 
TURBULENT  ENERGY 


ALL  Ax  POSSIBILITIES 


Figure  2.1  —  Flow  Logic  for  Subroutine  TURB 


25 


SSS-R-72-1255 


the  entire  grid  in  the  calculation  of  the  weights.  The  computer 
required  55  seconds  each  cycle  to  form  X2  and  J2  in  this  case. 

Because  of  the  great  expense  involved  in  carrying  out 
the  numerical  solution  of  the  turbulence  equations,  the  approach 
was  not  pursued  beyond  this  point.  It  is  expected,  however, 
that  some  effort  will  be  committed  during  the  next  year,  to  the 
modification  of  the  present  turbulence  scheme  to  speed  up  the 
computational  time  over  that  obtained  in  our  first  efforts. 

Only  then,  will  it  be  practical  to  attempt  calculation  of  the 
mountain  wave  problem  with  turbulence. 
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3.  LINEAR  STEADY  STATE  CODE  DEVELOPMENT 


The  calculation  of  the  vertical  flux  of  horizontal 

momentum  (wave  drag)  using  the  HAIFA  codes  as  described  pre- 
[121 

viously1  '  has  led  to  results  which  are  cyclic  with  time. 
These  results,  because  of  their  behavior,  are  not  well  under¬ 
stood  and  thus  have  not  been  parameterized.  In  order  to  more 
fully  understand  the  transient  results  and  to  aid  in  develop¬ 
ing  an  initial  parameterization  for  use  in  the  RAND  global 
circulation  model  (GCM) ,  a  linear  steady  state  numerical 
model  of  gravity  waves  was  developed.  It  is  important  to 
note  that  both  the  motion  of  the  flow  field  and  the  obstacle 
placed  in  the  stream  must  be  small  for  the  linear  code  to  be 
valid.  Some  of  the  problems  calculated  previously  using  the 
HAIFA  code  do  fall  in  this  category  and  it  was  felt  that  de¬ 
tailed  comparisons  could  be  made  with  reported  results  to 
aid  in  the  accomplishment  of  the  above  goals. 

The  code  description,  some  results  of  the  code  and 
a  parameterization  of  some  of  the  linear  steady  state  results 
are  described  in  this  chapter.  Comparisons  between  the  tran¬ 
sient  results  of  HAIFA  and  the  steady  state  results  are  also 
made  and  discussed  although  further  work  in  this  is  required 
before  the  transient  results  are  understood  well  enough  to 
parameterize  for  use  in  the  GCM. 


27 


SSS-R-72-1255 


3.1  LINEAR  STEADY  STATE  ANALYSIS 

Bretherton ^  and  Danielsen  and  Bleck, ^  among 
others,  have  published  linear  steady  state  analyses  on  the 
calculation  of  momentum  transport  by  gravity  waves.  After 
an  examination  of  both  procedures,  the  analysis  of  Brether- 
ton  was  chosen  as  the  basis  of  the  numerical  model  to  be 
developed  at  Systems,  Science  and  Software.  The  primary 
reason  for  the  selection  was  Bretherton  had  extended  his 
analysis  to  three  dimensions  and  arbitrary  topography  while 
the  Danielsen  and  Bleck  model  was  more  limited  in  scope. 

The  analysis  is  briefly  outlined  below.  The  reader  is  refer¬ 
red  to  the  referenced  papers  for  a  more  complete  discussion 
of  the  model. 

Small  perturbations  to  a  stratified  shear  flow  may 
be  analyzed  in  terms  of  a  single  Fourier  component  of  verti¬ 
cal  velocity 


w(x,z)  =  R  0(k,z)e^x  , 

where  k  is  the  wave  number  and  R  indicates  the  real  part 
of  the  complex  vertical  velocity  w(x,z).  Assuming  the  motion 
to  be  steady,  invisid  and  adiabatic,  making  the  Boussinesq 
approximation  and  neglecting  the  Earth's  rotation,  satisfies 
an  equation 


d2w  ,  ( 


dz : 


l2  (z)  -  k2|tf  =  0 


(3.1) 


r  7  1 

which  is  known  as  Scorer's  equation.  The  parameter  1  , 

sometimes  referred  to  as  Scorers  parameter  or  number,  is 
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. 2  _  N2 (z)  _  1  d2u(z)  n 

u2  (z)  u(z)  dz2 

r  81 

where  N  is  the  Brunt-Vaisala  frequency.  Eliassen  and  Pain1  1 
showed  that  the  Reynolds  stress  or  wave  drag  associated  with 
this  component  satisfies  the  equation 


fj(puw) 


d  (  p  /dw* 
dz  jzik \dz 


w 


dw  *  A ) 

“  d£  w*)j  = 


(3.3) 


where  p  is  the  mean  air  density,  and  w*  denotes  the  com¬ 
plex  conjugate  of  w  . 

The  equations  above  are  solved  subject  to  boundary 
conditions  at  the  lower  surface  and  at  the  top  of  the  atmos¬ 
phere.  The  lower  boundary  condition  is  based  upon  a  Fourier 
transformation  of  the  actual  topography  of  the  region  under 
study.  For  the  upper  boundary  condition,  Bretherton  divides 
the  problem  into  two  separate  parts  and  treats  each  part 
separately.  The  drag  associated  with  waves  which  are  assumed 
to  be  propagating  upwards  in  the  atosphere  above  the  arbitrary 
height  of  the  numerical  grid  (H)  are  calculated  separately 
from  those  associated  with  waves  which  are  trapped  with  the 
computational  atmospheric  grid. 

The  first  boundary  condition 


dw 

dz 


+  i/l2  -  k2  w  [sgr.  U] 


at  z  =  H 


(3.4) 


where  U  is  the  horizontal  free  stream  velocity,  is  applied 
when  k2  <  l2 (H)  and  implies  only  upward  propagating  waves 
above  Z=H.  The  second  boundary  condition, 
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^  =  -  Jk2  -  l2  w  at  Z  =  H  (3.5) 

dz 

is  applied  if  k2  >_  l2 (H)  and  implies  that,  should  the  basic 
state  of  the  atmosphere  be  imagined  to  continue  upwards  to 
infinity  with  l2  =  l2 (H)  ,  the  total  perturbation  energy  in 

this  wave  number  would  be  finite.  The  solution  to  this 
second  mode  (trapped  waves)  will  be  discussed  in  more  detail 
after  presenting  the  more  general  case  for  the  continuous 
drag  calculation. 

The  solution  of  Eq.  (3.3)  for  the  Reynolds  stress  is 
accomplished  in  the  following  manner: 

(1)  The  actual  surface  topography  is  Fourier  ana¬ 
lyzed  to  obtain  the  wave  number  dependence  of  the  lower 
boundary  condition. 

The  complex  Fourier  transform  is  defined  by 


fi  (k) 


-ikx 

e 


dx 


(3.6) 


where  h(x)  is  the  ground  level  height  of  the  position  x.  A 
spectrum  function  is  defined  as 

A(k)  =  h*h  (3.7) 

A 

where  X  is  the  maximum  grid  length  and  fi*  is  the  complex 
conjugate  of  h  . 
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(2)  The  equation  for  the  complex  velocity  w  ,  Eq. 
(3.1),  is  solved  subject  to  the  upper  boundary  condition, 

Eq.  (3.4) .  The  value  of  H  at  which  this  boundary  condition 
is  applied  is  somewhat  arbitrary  and  its  effect  is  investi¬ 
gated  in  test  problems,  the  results  of  which  are  presented 
later  in  this  section. 


From  the  results  of  the  integration,  a  parameter  F(k), 
independent  of  surface  topography,  may  be  computed  as 


F(k) 


1  jdw 
2T  jdz 


/v  .  ,  .  dw  *  /v  .  . 
W  ^  ~  dz~  w*z) 


w* (0)  w (0) 


(3.8) 


F (k) ,  which  is  independent  of  height,  has  the  property  that 
it  vanishes  if  l2 (H)  <  k2  .  To  obtain  the  entire  Reynolds 
stress,  we  add  the  effects  all  possible  wave  numbers,  with 
excitation  amplitudes  determined  from 


w(o)  =  ik  u£ 


(3.9) 


The  total  Reynolds  stress  from  Eq.  (3.3)  is  then 


puw  =  2 p  ( 0 )  U2(0)  /  k  A  (k)  F  (k)  dk 


(3.10) 


where  puw  =  Reynolds  stress,  U(0)  is  the  free  stream  hori¬ 
zontal  velocity  at  the  surface  (lower)  boundary,  and  p(0)  is 
the  fluid  density  at  the  lower  boundary.  It  may  happen  for 
a  range  of  values  kc  <  k  £  1(H),  that  F(k)  is  extremely 
close  to  zero  except  for  very  narrow  wave  number  intervals 
within  which  F(k)  becomes  large.  These  correspond  to 
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trapped  modes  with  a  very  slight  upwards  energy  leak.  Al¬ 
though  in  principle  F(k)  is  still  a  continuous  function 
of  k,  in  practice  it  would  be  very  inefficient  to  compute 
it.  numerically  from  Eq.  (3.10);  the  narrow  bands  where  F(k) 
differs  appreciably  from  zero  are  more  efficiently  treated 
analytically.  Consequently,  the  upper  limit  of  the  integral 
in  Eq.  (3.10)  is  reduced  to  kc  (defined  arbitrarily  as  the 
point  where  F(k)  approaches  zero)  and  the  trapped  wave  con¬ 
tributions  are  added  analytically. 

The  trapped  wave  analysis  outlined  in  Bretherton 
gives  for  the  Reynolds  stress 


puw  =  2  p ( 0 ) 


U2  (0) 


TT 
2 


AOc,) 


(3.11) 


where  the  subscript  i  is  the  ith  wave  number  greater  than 

k  but  less  than  £(H)  for  which  w(C)  =0  .  The  w.  is  the 
c 

solution  of  Eq.  (3.1)  with  the  boundary  condition  specified 
by  Eq .  (3.4). 

As  noted,  for  k  >  MH)  ,  F(k)  is  zero.  However, 
eigensolutions  satisfying  Eq.  (3.5)  with  w  =  0  at  the  ground 
may  exist;  in  this  case  the  lower  boundary  condition,  Eq. 

(3.9) ,  cannot  be  satisfied.  In  this  case,  F(k)  is  indeter¬ 
minate.  The  solutions  of  Eq.  (3.1)  subject  to  the  upper 
boundary  condition  specified  in  Eq.  (3.5)  are  found  and  the 
drag  is  calculated  using  Eq.  (3.11)  for  each  wave  number  for 
which  w(0)  =  0  . 
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3.1.1  Method  of  Numerical  Solution 

The  basic  scheme  used  in  the  solution  for  the  Reynolds 

stresses  associated  with  linear  steady  state  gravity  waves  is 

shown  in  Figure  3.1.1.  The  integrals  for  the  Reynolds  ft  cresses 

are  calculated  using  the  Gaussian  quadrature  method  while  the 

solution  for  the  complex  velocity  at  a  given  k  value  is 

found  using  a  Glauz-Adams  integration  ‘routine,  a  modification 

[9] 

of  the  Adams -Moulton  predictor-corrector  method. 

3.1.2  Test  Problems 

The  linear  steady  state  code  is  general  enough  to 

accept  any  temperature  stratification  and  horizontal  velocity 

profile  as  input.  However,  in  orler  to  compare  with  previous 

work,  the  problems  described  herein  have  constant  atmospheric 

lapse  rates  and  exponential  velocity  profiles  described  by 

U  eu  where  U  a:  .  UEXP  are  constants  for  each 

o  o 

problem.  In  this  way,  transient  results  of  problems  pre¬ 
viously  calculated  using  HAIFA  would  be  compared  with  the 

steady  state  results.  In  particular,  the  two  wave  and  single 

[1  2  ] 

wave  calculations  used  as  test  problems  '  which  were 
checked  using  the  analysis  of  Palm  and  Foldvik, could  be 
compared  with  the  results  of  this  new  code. 

A  matrix  of  problems  was  run  varying  the  atmospheric 
stability  and  the  exponential  velocity  coefficient,  using  a 
rectangular  mountain  625  meters  high  and  4.5  km  long.  H  was 
arbitrarily  selected  as  12.0  km,  roughly  corresponding  to  the 
grid  height  used  in  the  HAIFA  calculations.  The  effect  of 
this  height  on  the  results  is  discussed  later  in  this  sec¬ 
tion.  The  results  arising  from  each  calculation  are  (1)  the 
continuous  drag,  (2)  the  drag  associated  with  each  trapped 
wave,  (3)  the  trapped  wave  numbers  (wave  length) ,  and 
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Figure  3.1.1  —  Flow  chart  for  linear  steady  state  code 
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(4)  the  complex  velocity  perturbations  as  a  function  of 
height  in  the  atmosphere.  While  it  is  possible  to  compare 
the  wave  lengths  from  these  calculations  with  those  of  Palm 
and  Foldvik,  no  direct  comparison  is  available  of  the  drag 
values.  The  trapped  wavenumbers  are  tabulated  in  Table  3.1, 
as  a  function  of  the  exponential  velocity  coefficient  UEXP 
and  T' ,  the  atmospheric  lapse  rate  divided  by  T,  the  adiabatic 
lapse  of  the  atmosphere  (0.01  °C/m) . 

In  Figures  3.1.2  and  3.1.3,  the  real  part  of  the 
velocity  perturbation  as  a  function  of  altitude  is  shown  for 
two  problems.  The  first  is  the  two-wave  problem  resulting 
from  T'/r  =  0-5  and  UEXP  =  0.795  10-1*.  The  second  is  a 
problem  having  three  waves  defined  by  the  velocity  profile 
U  =  UqEXP  (1x10-1**Z)  and  an  atmospheric  lapse  rate  of 
0.002  °C/meter. 

The  drag  associated  with  each  wave  number  denoted  in 
Table  3.1  multiplied  by  the  wavenumber  is  presented  in  Table 
3.2  as  a  function  of  UEXP  and  T'.  The  general  trend  of  the 
results  indicates  larger  drag  values  as  (1)  the  lapse  rate 
decreases  compared  to  the  adiabatic  lapse  rate  and  (2)  the 
velocity  profile  approaches  a  uniform  value  (UEXP  =  0) .  The 
stronger  effect  appears  to  be  associated  with  the  lapse  rate. 
Continuous  drag  values  occur  only  at  low  values  of  the  velocity 
exponential  coefficient. 

The  two  wave  problem,  described  previously,  was  com¬ 
puted  with  a  grid  height  (H)  of  both  12  km  and  15  km.  The 
drag  associated  with  the  longer  wave  length  (K  =  0.1475)  de¬ 
creased  from  2.91  dynes/cm2  to  2.62  dynes/cm2.  The  wave 
number  itself  also  changed  by  approximately  2.5  percent 
while  the  shorter  wave  and  its  associated  drag  were  identical 
in  both  cases. 
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0 

u 

u 

TABLE  3.1 

WAVENUMBERS  (K)  OF  TRAPPED  WAVES 
AS  A  FUNCTION  OF 
LAPSE  RATE  AND  WIND  PROFILE 

L1 


u 

UEXP*104* 

(1/m) 

.  .  ** 
T'/r 

(lAm) 

K2  '(1/km) 

(1/km) 

L1 

0.5 

0.1 

1.421 

1.2037 

.1.0289 

1*5 

0.1 

1.0895 

0.6331 

0.2673 

1 1 

2.5 

0.1 

0.8146 

0.0748 

1 

i_j 

3.0 

0.1 

0.6787 

3.5 

0.1 

0.5354 

Li 

1.5 

0.2 

1.0074 

0.5624 

1.795 

0.2 

0.9246 

0.4142 

Lj 

1.0 

0.25 

1.1140 

0.7842 

0.5442 

2.0 

0.25 

0.8261 

0.2661 

3.0 

0.25 

0.5522 

| 

4.0 

0.25  . 

0.2265 

0.5 

0.40 

1.1450 

0.9581 

1 

1.0 

0.40 

0.9725 

0.6610 

— 1 

1.795 

0.40 

0.7447 

0.250 

, 

3.5 

0.40 

0.2392 

* 

— 1 

1.0 

J .  50 

0.8679 

0.5701 

1.795 

0.50 

0.6431 

0.1475 

I 

2.0 

0.50 

0.5875 

, 

1 _ J 

3.0 

0.50 

0.3688 

1.0 

0.75 

0.5475 

i 

2.0 

0.75 

0.2638 

, 

0.5 

0.90 

0.4138 

L 

1.5 

0.90 

0.1072 

Mountain  height  =  625  meters  Mountain  length  =  4.5  km 
*U  =  Uo  exp [UEXP*Z]  **r  „  0.01  oc/m 
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TABLE  3.2 

DRAG  (K)  AS  A  FUNCTION  OF  LAPSE  RATE, 
WIND  PROFILE,  AN  WAVE  NUMBER 


K1*DRAG(K1)  K2*DRAG(K2)  K3*DRAG(I<3) 
(dynes/cm2  x  km) 


0.0047 

1.407 

5.-354 


3.132 

0.312 


0.601 

1.097 


7.505 


1.795 


2.941 


2.709 

2.34 


0.736 

3.806 

5.826 

2.558 


1.879 

1.373 


1.795 


0.249 

1.213 

3.241 


0.578 

1.673 

0.982 


1.795 


1.425 

2.963 


1.399 

0.430 
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3. 1.2.1  Scaling  of  the  Results  —  The  drag  values  calculated 
from  the  linear  steady  state  mathematical  model  may  be  sepa¬ 
rated  into  a  continuous  value  and  a  value  associated  with  each 
trapped  wave  (or  angle  wave  number  k) .  Since  the  topography 
spectrum  enters  as  a  multiplicative  factor  in  the  drag  calcula¬ 
tion,  the  drag  values  in  the  tables  may  be  applied  to  any 
topography  by  multiplying  them  by  the  ratio  of  the  A(k)  for 
the  new  topography  data  to  A(k)  from  the  existing  runs.  For 
rectangular  obstacles  the  scaling  may  be  found  as  follows: 

A(k)  =  (3.12) 


where 


h  = 


1 

2  IT 


-ikx 

e 


dx 


(3.13) 


h(x)  corresponds  to  a  rectangular  mountain  for 
which  the  height,  h,  is  assumed  constant 
over  a  given  region  x.  to  x,  and  zero  else¬ 
where  in  the  grid. 


X  is  the  grid  length 


h  - h  r  -‘ 

h  -  27  / 

•'x. 


[coskx  -  i  sinkxjdx 


2irk 


[sinkx2  -  sinkx1  +  i(coskx2  -  coskx^  ]  (3.14) 


h*  = 


2?rk 


[sinkx2  -  sinkx1  -  i(coskx2  -  coskx^ ]  (3.15) 


A(k)  =  x2iTk j  t(sinkx2  -  sinkx^) 2  +  (coskx2  -  coskx^)2] 

(3.16) 
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Thus,  we  note  that  for  rectangular  mountain,  the  drag  for 
any  single  wave  number  is  proportional  to  the  height  squared 
and  to  the  mountain  length  in  a  more  complicated  function  of 
the  square  of  the  sum  of  sines  and  cosines.  The  formula  can 
be  further  simplified  by  noting  that  x2  -  x1  =  L  where  L 
is  the  mountain  length.  Making  this  substitution,  A(k)  = 
(h2/Xirk2)  (1-coskL)  . 

The  results  of  the  calculations  have  been  parameterized 
so  that  (1)  the  primary  and  secondary  wave  numbers  may  be 
found  as  a  function  of  UE7P  and  T'/r  and  (2)  the  drag  associ¬ 
ated  with  each  wave  may  be  calculated  from  the  same  two  param¬ 
eters.  Figures  3.1.4  and  3.1.5  present  the  wave  number  cal¬ 
culations  and  contours  selected  to  represent  the  results.  The 
equation  for  k^  (the  primary  wave  number)  was  found  to  be 

k1  =  — 0 . 32  8  [  (UEXPxlO1*)  +  4(T'/r)]  +  1.78  (1/km)  (3.17) 

and  k2  (the  secondary  wave  number)  as 

k2  =  -0 . 431  [  (UEXPx.10 4 )  +  2(T'/r)l  +  1.48  (1/km)  .  (3.18) 

While  some  of  the  cases  calculated  showed  the  presence  of  a 
third  wave,  an  insufficient  number  of  points  were  available 
to  parameterize  the  results.  It  is  important  to  point  out 
that  the  above  equations  may  be  improved  upon  as  more  data 
from  the  calculations  become  available.  Also,  it  was  deter¬ 
mined  that  in  general  a  negative  value  of  k^  or  k2  indicates 
no  wave  is  present  in  that  region. 

The  k^ (Drag (k^) /A(k^) )  and  k2 (Drag (k2 ) /A(k2) )  are 
shown  in  Figures  3.1.6  and  3.1.7.  These  values  were  also 
parameterized  as  functions  of  UEXP  and  T'/r  and  can  be  repre¬ 
sented  by  the  equations: 
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k1|Drag(k1)/A(k1)J  =  ^j-  6.45  +  (UEXPxlO1*) 

±  ^-HxUEXPxlO4  -  16.8x11+  54.2j1/2j  (3.19) 


k2|Drag(k2)/A(k2)J  =  -  9.44xl02  Jl  +  ‘7.26xi02 

-  2 . 36x10 2  j UEXPxlO  4  -  0 . 85  J 


(3.20) 


The  actual  values  of  k(Drag(k)/A(k) )  used  for  the 
parameterization  are  included  in  Table  3.3  As  previously  noted, 
the  drag  can  bo  scaled  for  any  topography  for  which  A(k)  can 
be  calculated.  Further  efforts  will  be  made  to  present  tht 
results  in  the  manner  most  useful  for  incorporation  into  the 
GCM.  In  addition,  the  effects  of  using  wind  profiles  other 
than  exponential  in  shape  and  the  effects  of  non-constant  at¬ 
mospheric  lapse  rate  will  be  investigated. 


3.2  COMPARISON  OF  NON-LINEAR  TIME  DEPENDENT  AND  LINEAR 

STEADY  STATE  SOLUTIONS 

For  two  problems,  a  comparison  between  the  linear 
steady  state  code  and  the  time  dependent  HAIFA  code  can  be 
made.  These  were  the  single  wave  and  two  wave  problems  re¬ 
ported  in  Reference  1  and  which  also  have  been  compared  with 
the  results  of  Palm  and  Foldvik's  analysis.  The  problem  in¬ 
put  is  characterized  in  the  following  table. 

Figure  3.2.1  presents  the  results  of  the  drag  as  cal¬ 
culated  using  the  two  codes.  It  is  seen  that  the  HAIFA  re¬ 
sults  would  have  to  be  compared  to  much  later  time  in  order  to 
determine  or  justify  an  extrapolation  of  the  results  to  the 
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TABLE  3.3 

CALCULATED  VALUES  OF  k (Drag (k) /A (k) ) 
USED  IN  THE  PARAMETERIZATION 


UEXP 

T'/ 

k2  (Drag/A (k2))  . 

k^  (Drag/A (k^) ) 

1.5 

0.2 

378.35 

1335.7 

0.5 

0.25 

475.9 

707.96 

1.5 

0.35 

212.3 

885.5 

1.0 

0.4 

295.6 

697.1 

2.0 

0.6 

267.9 

2.0 

0.7 

124.6 

1.0 

0.65 

84.76 

250.6 

3.5 

0.5 

41.21 

3.0 

0.6 

58.2 

3.0 

0.1 

1299.0 

2.5 

0.3 

848.4 

2.5 

0.45 

457.4 

1.795 

0.2 

249.7 

1336.8 

1.795 

0.4 

86.1 

731.6 

1.795 

0.6 

1.7 

296.3 

1.795 

0.9 

46.2 

1.0 

0.25 

206.0 

483.0 

2.0 

0.25 

122.6 

1139.0 

3.0 

0.25 

796.0 

4.0 

•  0.25 

220.5 

1.0 

0.75 

129.5 

2.0 

0.75 

• 

70.6 

3.5 

0.1 

990.1 

2.5 

0.1 

25.0 

1534.0 

1.5 

0.1 

98.0 

515.0 

0.5 

0.1 

477.0 

.  678.0 

3.5 

0.4 

185.0 

0.5 

0.4 

309.0 

465.0 

0.5 

0.7 

131.0 

2.5 

0.6 

176.0 

1.5 

0.9 

7.0 

0.5 

0.9 

17.7 
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PROBLEM 

ATMOSPHERIC 
LAPSE  RATE 

OBSTACLE  SIZE 

HORIZONTAL  VELOCITY  PROFILE 

Single  Have 

1/2  r111 

635  meters  high 

U  - 

9.43 

exp  (3  «  10"'  z)  121 

by 

4.5  km  long 

U  - 

10.0 

exp  (3  x  io"'  z)  131 

‘  ^  V 

Two  Have 

1/2  r 

625  meters  high 

U  - 

9.94 

exp  (1.795  f  10"'  z) 

[2] 

by 

4.5  km  long 

U  - 

10.0 

exp  (1.795  x  10“'  z) 

[3] 

^  r  «  adiabatic  lapse  rate  in  the  atmosphere 


11  HAIFA  code 

^  Linear  steady  state  code 


steady  state  value.  This  agrees  with  the  previous  interpre¬ 
tation  of  the  transient  results  in  that  the  Reynolds  stress 
was  far  from  achieving  a  constant  value  as  a  function  of  height, 
a  requirement  for  the  steady  state  case. 

3.3  COMPARISON  OF  RESULTS  OF  THE  LINEAR  STEADY  STATE 

ANALYSIS  AND  A  GCM  PARAMETERIZATION 

The  parameterization  of  the  drag  currently  used  in 
the  Rand  Climate  Dynamics  version  of  the  UCLA  Global  Circu¬ 
lation  Model ^  ^  relates  the  surface  stress  (t  )  to  the  total 

s 

surface  velocity,  U  ,  and  a  drag  coefficient,  Cn,  by  the  rela- 

S  l) 

tion 


Ts  =  -P  CD  |Us  +  G|  Us  (3.21) 

where  p  is  the  air  density,  U  is  estimated  by  a  linear 

b 

extrapolation  of  U  from  the  two  known  levels  with  respect 
to  the  vertical  pressure  coordinate  and  G  is  a  gustiness 
factor  taken  to  be  2.0  m/sec.  The  drag  coefficient  over  land 
is  taken  to  be 
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CD  =  0.002  +  0.006  Zg/5000  (3.22) 

where  Zg ,  the  average  altitude  of  the  surface,  is  in  meters. 

The  mountain  range  used  in  the  linear  steady  state 
calculations  was  625  meters  high.  The  drag  coefficient,  CD, 
in  these  circumstances  (Eq.  (3.  ) )  is  0.00275  producing  a 

constant  drag  of  3.30  dynes/cm2  for  a  -surface  wind  velocity 
of  10.0  meters/sec. 

The  results  calculated  from  the  linear  steady  state 
model  indicate  differences  that  should  be  included  in  the 
parameterization  in  at  least  two  different  areas.  First,  the 
results  indicate  a  dependence  of  the  drag  on  the  horizontal 
velocity  profile.  Secondly,  the  variation  of  the  mountain 
width  and  the  atmospheric  lapse  rate  calculated  produced  a 
variation  in  drag  values  which  ranged  from  zero  to  values 
greater  than  10.0  dynes/cm2. 

One  concludes  from  the  very  limited  studies  completed 
that  an  average  drag  based  on  the  surface  roughness  alone  is 
not  adequate  to  account  for  orographic  effects  of  the  vertical 
flux  of  horizontal  momentum.  The  effects  of  lapse  rate 
mountain  widths,  and  horizontal  velocity  profiles  should  also 
be  addressed.  Further  calculations  in  both  two  and  three 
dimensions  will  be  carried  out  in  order  to  determine  in  a 
more  definitive  manner  how  the  wave  drag  is  affected  by  these 
variables . 
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4.  RADIATION  IN  THE  EARTH'S  ATMOSPHERE: 

ATRAD  THEORY  * 

Since  the  previous  semiannual  report  on  this  con- 
[2] 

tract,  the  atmospheric  radiation  code,  called  ATRAD,  has 
become  operational.  It  has  been  compared  against  several 
known  radiative  transfer  solutions  and  the  agreement  has  been 
excellent.  Two  comparison  calculations  have  been  made  be¬ 
tween  ATRAD  and  the  radiation  subroutine  of  the  Mintz-Arakawa 
code.  The  results  of  these  comparisons  and  a  description  of 
ATRAD  were  presented  at  the  American  Meteorological'  Society ’ s 

Conference  on  Atmospheric  Radiation  at  Fort  Collins,  Colorado, 

[121 

August  7-9,  1972,  1  (This  paper  is  included  as  Appendix  C.) 

Since  the  code  ATRAD  has  reached  operational  status 
and  fewer  changes  will  be  made  in  the  future,  it  is  useful  now 
to  summarize  its  theoretical  basis.  This  is  done  in  the  pre¬ 
sent  chapter,  where  elements  from  both  of  the  prior  semi- 
[1  21 

annuals  '  and  code  improvements  which  have  not  previously 
been  reported  are  brought  together.  Chapter  5  contains  a 
description  of  the  Mintz-Arakawa  radiation  subroutine  and  of 
the  method  of  obtaining  data  from  it  to  form  appropriate  input 
data  for  ATRAD.  The  two  comparisons  between  ATRAD  and  Mintz- 
Arakawa  are  also  presented  and  discussed  in  Chapter  5. 

Finally,  future  plans  for  code  development  and  code  applica¬ 
tions  are  presented  in  chapter  6 . 


51 


SSS-R-72-1255 


4.1  ASSUMPTIONS 

The  fundamental  assumptions  involved  in  the  formula¬ 
tion  of  the  ATRAD  equations  are: 

(1)  plane-parallel  atmosphere  (rather  than 
spherical)  and  horizontal  surface, 

(2)  horizontal  homogeneity, 

(3)  local  thermodynamic  equilibrium,  so  that 
true  emission  is  described  by  the  Planck 
function, 

(4)  unpolarized  radiation  field  , 

(5)  non-re tractive  air  (rays  are  not  curved) , 

(6)  spherical  aerosol  scatterers  of  uniform 
composition , 

(7)  point  source  sun  (no  angular  width) ,  and 

(8)  multiplying  together  the  transmissions  of 
individual  molecular  constituents  yields 
the  total  transmission. 

Assumptions  (1)  and  (5)  are  only  violated  for  large- 
zenith-angle  rays,  and  such  rays  make  only  a  very  small  con¬ 
tribution  to  vertical  fluxes. 

Assumption  (2)  probably  represents  the  most  serious 
approximation  in  the  code,  one  whose  impact  has  been  largely 
ignored  in  the  literature.  The  effect  of  the  most  important 
form  of  horizontal  inhomogeneity,  partial  cloudiness,  is  cus¬ 
tomarily  accounted  for  by  taking  a  weighted  sum  of  clear-air 
and  totally-cloudy  fluxes, 

F  =  aF  ,  ,  +  (l-a)F  , 

cloud  clear 
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where  a  is  the  fractional  cloudiness.  The  only  virtue  of 
this  scheme  is  that  it  reduces  properly  in  the  limits  a  =  0 
and  a  =  1  and,  therefore,  cannot  err  too  greatly  for 
0  <  a  <  1  .  A  more  reasonable  scheme,  in  the  context  of 
ATRAD's  numerical  algorithm, ^  would  be  to  take  the  weighted 
sum  of  the  reflection  and  transmission  matrices  of  a  partially 
cloudy  zone.  This  difficulty  could  be  resolved  with  a  three- 
dimensional  Monte  Carlo  code  such  ac  that  of  Cox  and  McKee, 
which  should  help  to  ascertain  the  relative  importance  of 
horizontal  inhomogeneity  and  whether  it  is  really  significant 
for  weather  and  climate  radiation  subroutines. 

Assumption  (3) ,  which  requires  sufficiently  high 
pressures  that  collisions  maintain  thermodynamic  equilibrium 
in  excited  states,  is  valid  below  about  70  km. 

Assumption  (4)  has  been  examined  by  several  investi¬ 
gators,  for  example  Howell  and  Jacobowitz ,  who  find  that, 
especially  in  the  presence  of  aerosols,  the  flux  errors  re¬ 
sulting  from  a  polarization-independent  treatment  are  generally 
less  than  one  percent  in  the  visible  spectrum.  Averaged  over 
the  whole  spectrum,  these  errors  are  considerably  smaller. 

.  r  1 4 1 

A  standard  Mie  theory  treatment  of  aerosol  scatter¬ 
ing  is  made  possible  by  assumption  (6) .  Uniform  composition 
is  not  essential,  since  layered  spherical  particles  can  be 

handled  with  some  increase  in  computational  complexity  (see 
[15] 

Kerker  ).  However,  there  is  practically  no  experimental 
data  on  layered  (water-sheathed,  etc.)  aerosol  particles. 
Neither  is  sphericity  essential  since  Mie  theory  has  been  ex¬ 
tended  to  ellipsoids  and  circular  cylinders.  However,  solid 
aerosol  particles  (with  the  exception  of  cirrus  ice  needles) 
are  no  more  accurately  representable  as  ellipsoids  or 

* 

S.  Cox,  Colorado  State  University,  Fort  Collins,  Colorado, 
private  communication. 


53 


SSS-R-72-1255 


cylinders  than  as  spheres.  And  liquid  aerosol  particles  are 
very  nearly  spherical,  although  the  larger  ones  will  tend  to 
be  slightly  flattened  by  the  air  flow  around  them.  Further¬ 
more,  as  Hodkinson  ^  ^  argues,  randomly  oriented  irregular 

particles  produce  the  same  forward  diffraction,  the  same 
scattering  due  to  external  reflection,  and  the  same  scatter¬ 
ing  due  to  the  first  refraction  as  equivalent  spheres  (spheres 
having  the  same  distribution  of  projected  area) .  For  moder¬ 
ately  absorbing  irregular  particles,  this  x^ould  lead  to  a 
complete  equivalence  of  scattering  patterns.  Measurements 
of  Hodkinson, 116 ]  Ellison,  [17]  and  Holland  and  Gagne [18] 
tend  to  confirm  the  usefulness  of  Mie  theory  (with  equivalent 
spheres)  for  predicting  the  scattering  from  irregular 
particles  into  the  forward  hemisphere,  which  contains  almost 
all  the  scattered  photons. 

The  one  aerosol  which,  by  virtue  of  its  importance 
to  the  global  radiative  energy  budget, 119 ' 20 ^  should  probably 

not  approximated  as  spherical  is  cirrus  cloud. 

[21] 

Jacobowitz  has  computed  the  phase  function  of  randomly 
oriented  hexagonal  ice  cylinders,  typical  of  cirrus,  and  has 
found  that  as  much  energy  is  scattered,  into  the  halo,  around 
22°,  as  into  the  main  diffraction  peak.  This  effect  is  en¬ 
tirely  missed  by  an  equivalent  sphere  model.  In  the  future, 
we  intend  to  include  cirrus  in  ATRAD  more  completely. 

The  error  incurred  by  assumption  (7) ,  ignoring  the 
angular  width  of  the  sun,  has  been  estimated  and  found  to  be 
negligible  as  far  as  vertical  fluxes  are  concerned. 

Assumption  (8),  that  individual  transmissions  are 
multiplicative,  is  firmly  grounded  in  experiment,  and 

lias  virtually  no  cetectable  experimental  error  associated 
with  it. 
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4.2  BASIC  RADIATIVE  TRANSFER  EQUATIONS 


The  radiative  transfer  problem  in  the  Earth's  atmos¬ 
phere  reduces  to  the  solution  of  the  seemingly  simple  equa- 

..  [24] 

tion 


"v  .  .  .  K  r 

ds  Jv  Kv  v' 


(4.1) 


which  states  that  the  monochromatic  (frequency  v)  radiant  in¬ 
tensity  I  ,  in  traversing  the  element  of  length  ds  ,  will  be 
augmented  by  sources  in  the  amount  ds  and  diminished  by 
extinction  in  the  amount  I  ds  .  In  general,  the  radiant 

intensity  and  J  ,  the  source  function,  depend  on  both  a 
spatial  coordinate  r  and  an  angular  coordinate  (direction) 
fi  at  the  point  r  ,  as  well  as  upon  the  frequency  v  .  The 
time  dependence  of  these  quantities  is  ignored  because  the 
radiative  state  of  the  atmosphere  is  established,  for  all 
practical  purposes,  instantaneously,  is  the  extinction 

coefficient,  which  describes  the  relative  depletion  in  the 
intensity  of  the  beam,  dl^/I^  ,  upon  traversing  the  element 
of  distance  ds  .  is  in  general  the  Siam  of  an  absorption 

part  ctnd  a  scattering  part.  describes  the  additions  made 

to  the  beam  intensity  along  ds  by  thermal  (Planckian)  emis¬ 
sion  and  by  scattering. 

In  ATRAD's  plane-parallel  geometry,  which  is  illustrated 

in  Figure  4.1,  r  is  replaced  by  z,  which  is  measured  from 

the  top  of  the  atmosphere,  and  is  represented  as  usual  by 

the  spherical  coordinates  y  =  cos9  and  <j>  .  Using  Kirch- 

ftoff's  Law  to  obtain  the  source  function  J  ,  the  radiative 

V  ros  1 

transfer  equation  (4.1)  becomes  in  this  geometry 
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77V 


(surface) 


Figure  4.1  -  ATRAD  coordinate  system 


av ( z)Bv (T)  + 


-  Kv<z> 


5^^-  f  P,(z,M')  Ijz^'jdn' 

4tt  /  V  v 

<4  TT 


=  a'  +  8  =  volume  extinction  coefficient, 

v  v 

=  volume  absorption  coefficient,  ,  corrected 
for  stimulated  emission  (a  correction  usually 
ignored  in  the  atmospheric  radiation  literautre) , 


L  -hv/kT\ 
=  %  1  "  e 


r  =  temperature, 


=  volume  scattering  coefficient, 


=  Planck  function 
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and  is  the  phase  function,  defined  so  that 


P  yj  (  Z  ,  ^  ^  ) 


dfi 

4tt 


is  the  probability  that  a  photon  entering  a  volume  element 
around  z  from  direction  will  be  scattered  into  the 

cone  dfl  of  directions  around  ^  .  Since  absorption  is  ex¬ 
plicitly  represented  in  Eq.  (4,2),  the  above  probabilities 
must  sum  to  one 


(4.3) 


In  what  follows ,  it  is  convenient  to  break  the  absorp¬ 
tion  coefficient  into  a  "line"  part  and  a  "continuum"  part, 


„ i  _  line  ,  cont 
a  =  a  +  a 

v  v  v 


(4.4) 


cont 

av  vanes  slowly"  with  frequency  and  includes,  for  example, 

the  water  vapor  8-12y  continuum  and  the  aerosol  absorption; 
line 

«v  ls  the  rapidly-varying  part,  due  to  absorption  lines  of 
atmospheric  constituents .  Since  a  continuum  is  often  due  to 
the  wings  of  distant  lines ,  the  separation  is  to  a  certain  ex¬ 
tent  arbitrary.  Therefore,  we  shall  adopt  the  convention 
that,  with  the  exception  of  aerosol  absorption,  contjnua  exist 
only  in  the  gaps  between  absorption  bands. 

Because  a  volume  element  in  the  atmosphere  exhibits  the 
same  scattering  pattern  no  matter  what  the  direction  of  incidence 
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,  the  phase  function  depends  only  on  the  scattering  angle 

0  between  ^  and  , 
s  - — 


P 

V 


where 


y  _  =  cos0  =  yy'  +  /l-y2  A-y'2  cos  ( 4>— <J> ' ) 

s  s 


In  order  to  azimuthally-average  the  radiation  transport  equa¬ 
tion  (see  below),  it  is  necessary  to  azimuthally-average  P^: 


Pv(z,y,y') 


_  1_ 

'  2  IT 


1 

2ir 


Pv(z,ys)d()) 


Pv(z,yy'  + 


A-y  2  A-y '  2  cos<j>)d(J> 


(4.5) 


The  second  step  follows  from  the  periodicity  of  the  cosine. 
As  a  result,  P^  does  not  depend  on  <j>'  . 

For  the  computation  of  vertical  fluxes  and  radiative 
heating  rates,  the  retention  of  the  (^-dependence  of  I  is 
unnecessary.  It  can  be  eliminated  by  ^-averaging  Eq.  (4.2) 

i  r2v 

(operating  on  both  sides  with  /  d0)  to  obtain 
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31 

y  3^  =  3v  + 


f1 

2  J- 1 


p  (z,u,y')  I  (z,y')dy*  -  k  I 


V  V 


(4.6) 


where 


Iv(z,y)  = 


=  ~  f 

27r  Jo 


2t r 


I  (z,y,<f>)d<(>  . 


Knowing  Iv  ,  the  vertical  flux  F  is  obtained  from 


Fv(z)  = 


-  f  y  I.,  (z  ,^)  dfl  =  2 TT  f 

j  4*  V  J- 1 


y  Iv(z,y)dy 


(4.7) 


From  Fv  the  total  (spectrally-integrated)  flux  F  can  be 
computed: 


F (z)  =  f  max  F  (z)dv 

A>  .  v 
mm 

where  ^ vmin ' Vmax^  includes  all  spectral  intervals  in  which 
Fv  is  n on-negligible .  From  F  ,  the  radiative  heating  rate 
at  level  z  can  be  calculated  from 

9T  _  1  dF 

"  pc  dz  (4.8) 

P 


where 


P  -  p(z)  =  air  density 

p  -  specific  heat  of  air  at  constant  pressure  . 
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Using  the  hydrostatic  approximation,  this  becomes 

9T  _  dF 

Tt  “  cp  dP  * 

Rather  than  compute  point-values  of  heating,  ATRAD  computes 
the  average  heating  of  a  zone  (z^z^  ^)  in  terms  of  the  bound¬ 
ary  fluxes: 


(4.9) 


(the  slight  vercical  variation  of  g  has  been  ignored) . 

The  radiative  transfer  equation,  Eq.  (4.6),  is  a 
phenomenological  equation  requiring  input  data  of  three  dif¬ 
ferent  types : 


(a)  atmospheric  structure  data 

(b)  boundary  condition  data 

(c)  absorption  data  for  relevant  atmospheric 
constituents . 


The  next  three  sections  describe  ATRAD 1 s  needs  in  each  of  these 
areas . 
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4.3  SPECIFICATIONS  OF  ATMOSPHERIC  STRUCTURE 

The  atmospheric  structure  data  needed  to  solve  Eq. 
(4.6),  with  the  use  of  Mie  scattering  theory  for  aerosols, 
consists  of  the  pressure  p,  the  temperature  T,  the  water 
vapor  density  p  ,  the  ozone  density  p__,  and  the  aerosol 
(including  cloudf  number  density,  size  distribution,  and  compo¬ 
sition.  These  must  all  be  given  as  a  'function  of  height. 

The  atmospheric  absorbers  other  than  water  vapor  and  ozone 
are  assumed  to  be  uniformly  mixed. 

The  details  as  to  how  the  structure  data  are  sup¬ 
plied  to  ATRAD  are  given  in  Appendix  A.  A  variety  of  analytic 
aerosol  size  distributions  are  available  as  options,  or  the 
user  may  supply  his  own  size  distributions  in  card  form, 
similarly ,  a  variety  of  aerosol  refractive  indices  are'  avail¬ 
able  as  options,  e.g.  water,  ice,  sea-salt,  dust,  etc. 

4 . 4  BOUNDARY  CONDITIONS 

The  boundary  condition  at  the  top  of  the  atmosphere 
is 


Iv(0,U,<f>)  =  Sv  6  (ft  - 


0  <  y  <  i 

for 

0  <  <J>  <  2tt 


where  Sv  is  the 
the  Earth  and  ^ 

o 

coordinate  system 


solar  flux  at  frequency  v  at  the  position  of 
is  the  direction  of  the  solar  beam  in  ATRAD' s 
(Figure  4.1)  .  Taking  the  4>-average, 


Iv(0,u)  -  Ml*  -  V 


0  <  y  <  1 


(4.10) 
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yo,  the  cosine  of  the  solar  zenith  angle,  is  calculated  follow¬ 
ing  the  prescription  of  Woolf f26]  from  user-supplied  values  of 
latitude,  longitude,  day  of  the  year,  and  Greenwich  time.  The 
best  measurements  of  Sv  are  those  of  Thekaekara,  et.al.,l27] 
which  have  been  included  in  the  code.  (Actually,  the  values  of 


dv 


are  given  by  Thekaekara,  and  the  code  obtains  from  these  the 
values 


Sav  Av  Xv 

for  use  in  each  frequency  group  Av.)  SAv  is  adjusted  accord¬ 
ing  to  the  earth-sun  distance,  which  is  calculated  from  the 
day  of  the  year. 

The  boundary  condition  at  the  surface  requires,  for 
its  complete  specification,  the  ground  temperature  T  and 
the  bidirectional  reflectivity  pv".  if  a  beam  of  intensity 

Iv,inc  is  incident  on  a  surface  from  direction  ifj.  and  the 
resultant  reflected  intensity  at  angle  <5  is  I1  (ft  jtl.), 
then  pv"  is  defined  such  that[28]  (see  Figure  4. 2)^  r  1 


TTI 


v,ref(fir;^i}  "  pv  {^i'\);Ev,inc  coseidfii  (4.11) 


By  this  definition,  is  of  finite  order,  barring  specular 

reflection,  for  I^ref  is  of  differential  order  with  respect 

t0  , inc ’  In  case  there  is  a  specular  component  of  reflec¬ 
tion,  we  separate  p"  into  a  diffuse  part  p"  and  a 

V  |  U 
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specular  part 


P 


ii 

v,s 


/ 


P 


H 

\> 


P 


v,d 


+ 


where 


p"  (h .  ,-fo  )  =  —2L_ 

Mv,  s  '  i  r  cos 6 


pv  s^i)6  ^yi“yr)<S  (<^i+Tr“<f)r) 


(4.12) 


Some  authors,  for  example  Siegel  and  Howell^28-'  ,  do  not 
include  the  factor  of  it  in  Eq.  (4.11).  We  find  it  convenient 
because  when  the  reflection  is  diffuse,  reduces  to  an 

albedo  (flux  ratio) ,  the  customary  measure  of  surface  reflecti¬ 
vity. 


Figure  4.2.  Geometry  of  reflection. 
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When  radiation  is  incident  from  all  1h,  the  reflected 
intensity  along  tj  is,  by  Eq.  (4.11), 


I  (fl  )  =  ± 

v,ref  r 


i  rzv  r 1 

=  rf  / 

o  r 


d^i 


v 


(ft.  )]J. 

,inc  i  r 


Adding  on  the  thermal  emission  from  the-  ground  and  phrasing  the 
whole  result  in  terms  of  the  ATRAD  coordinate  system  (Figure  4.1), 
the  surface  boundary  condition  becomes 


Iv ( z  o • P • $ ) 


£;(|y|)Bv(Tg)  +  p^s(|y 

|)IV  (z0 • ! y  1  # 4>-tt ) 

i  r2TT 

+  ?J  d*’ 

J  o  * 

1  jd|J,pv,a!  I*1 

1 , y  '  # '  )IV  (zfl  ,y 

(4.13) 
<j>'  )y ' 


for 


-  1  £  y  <  0 
0  <  (j)  <  2tt 


where  zg  is  the  surface  level,  and 


c  * 

ev' 


the  directional 


emissivity,  has  been  written  as  a  function  of  angle  of  emis¬ 
sion.  The  absolute  value  signs  on  y  in  e’  o’  .  and 

v '  ^ v,s ' 

Pv,d  are  merely  because  these  surface  properties  are  customarily 
measured  or  calculated  with  the  convention  0  <_  0^,  6  <  90°, 

whereas  ATRAD' s  coordinate  svstem  takes  0  £  0^  <  90°  and 
0r  <  180°. 


90°  < 


It  should  be  noted  that,  in  Eq.  (4.13),  <j>  is  entirely 

missing  from  the  argument  list  of  p'  and  occurs  in  the  com- 

V  f  s 

bination  <j>-<f>'  in  the  argument  list  of  p"  .  In  both  cases, 

v  f  d 

this  is  justified  because  the  surface  is  isotropic,  that  is. 
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its  reflecting  properties  do  not  change  as  we  rotate  the  co¬ 
ordinate  system  in  Figure  4.2  about  the  vertical  axis. 

Taking  tie  ^-average  of  Eq.  (4.13)  leads  to  the  final 
form  of  the  surface  boundary  condition: 


Vzo'^  =  ^CybVTg)  +  p;f8(lnl>vv  lyD 


i, 


(4.14) 


+  2J  dy  1  Pv/d{  b  l»y’  )Iv(z0,y')y' 
0 


for  -  1  <  y  <  0  ,  where 


-2TT 

PC,d(lJ'y,)  =271 

0 

27T 

Pv,d  {y,y' 

o 

(the  second  step  follows  from  the  2ir-periodicity  of  p"  ,  in 

V  ^  u 

its  third  argument) . 

There  are  three  simplifications  of  the  surface  boundary 
condition  (4.14)  which  are  used  in  almost  all  of  the  other  ex¬ 
tant  atmospheric  radiation  codes: 

l'a)  P;,s  =  pv,d  =  °‘  This  is  used 

practically  without  exception  by  codes 

dealing  with  the  infrared  (IR)  spectrum, 
ignoring  the  very  real  frequency-^and 
angle-variations  of  for  natural  sur¬ 

faces.  Furthermore,  except  for  water, 
snow,  and  ice,  the  average  IR  emissivity 
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of  natural  surfaces  is  known  to  deviate 

[221 

considerably  from  unity*  Nevertheless, 

the  mathematical  consequences  of  using  a 
non-unit  emissivity  (implyinq  p^  ^  ^  0) 
are  disastrous  for  codes  using  a  simple 
quadrature  of  the  integral  transport  equa¬ 
tion,  for  then  they  can  no  longer  frequency- 
average  the  boundary  term: 

pv,d  =  °'  pv,s  given  by  Fresnel  formulae 
for  reflection  at  a  plane  interface,  v  so 
large  that  B^(Tg)  is  negligible.  This  is 
used  primarily  by  codes  dealing  with  single 
frequencies  in  the  solar  spectrum  with 
water  as  the  underlying  surface. 

p;,s  "  °'  pv,d  =  constant  (independent  of 

angle  of  incidence  and  angle  of  reflection) , 

v  so  large  that  Bv(Tg)  is  negligible. 

This  is  used  primarily  by  codes  dealing  with 

single  frequencies  in  the  solar  spectrum 

with  a  solid  underlying  surface.  p"  ,  is 

V  f  Cl 

the  albedo,  and  the  reflection  is  assumed 
diffuse. 

We  believe  that  it  is  important  to  include  a  fully  general 
boundary  condition,  Eq.  (4.14),  as  ATRAD  does.  Research  to 
date  suggests  that  such  factors  as  the  angular — and  frequency — 
dependence  of  or  p^  cannot  be  ignored. [28 ' 30 ' 31 ' 32 » 33 ) 

While  data  with  respect  to  these  dependences  are  limited,  they 
exist  in  sufficient  quantities  to  warrant  comparisons  with 
calculations  using  frequency  and/or  angu’ ar  averages.  Such 
studies  will  be  performed  with  ATRAD,  and  should  indicate 
what  level  of  detail  is  necessary  in  the  boundary  condition. 
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Two  theoretical  models  of  rough-surface  reflectivity 
seem  sufficiently  realistic  to  warrant  inclusion  in  ATRAD. 

The  first  is  that  for  a  rough  water  surface,  in  which  it  is 
assumed  that  the  surface  consists  of  planar  elements  each  of 
which  is  a  specular  reflector.  The  slope  distribution  of  the 
Planar  elements  as  a  function  of  average  surface  wind  speed 
can  be  taken  from  the  experimental  work  of  Cox  and  Munk. [29] 
Even  with  this  model,  the  inclusion  of  shadowing  and  multiple 
reflections  among  planar  elements  is  difficult.  Therefore 
these  effects  can  probably  only  be  included  approximately. 
Experimental  measurements  of  sea  surface  bidirectional  re¬ 
activities  should  prove  useful  .n  guiding  the  ima_ 

tions. 


The  second  theoretical  model  of  rough-surface  reflecti¬ 
vity  is  for  land  areas,  and  assumes  that  the  surface  consists 
or  planar  elements  each  of  which  is  a  diffuse  reflector  (of  given 
abedo).  This  model  is  due  to  HcStraviok.  1311  The  slope  distri¬ 
bution  of  the  planar  elements  is  presumed  Gaussian.  As  with  *he 
first  model,  shadowing  and  multiple  reflection  can  probably  be 
included  only  approximately. 


In  practice,  the  full  bidirectional  reflectivity  of  the 
earth  s  surface  is  not  available  experimentally,  nor  is  it 
likely  to  be  available  in  the  near  future,  as  a  function  of 
latitude  and  longitude.  Therefore  the  two  theoretical  models 
described  above,  particularly  the  water-surface  model  applying 
to  so  much  of  the  earth's  surface,  assume  special  importance 
in  filling  this  void.  However,  ATRAD  also  makes  provision  to 
use  certain  degraded  forms  of  the  bidirectional  reflectivity 
and  directional  emissivity  which  are  more  readily  available 
from  experiment.  These  are  (refer  to  Figure  4.2  for  notation): 


(a)  the  uirectional-hemispherical  reflectivity. 
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•+  reflected  flux 

(^i )  — 

v  ■*-  incident  flux  from  u. 


f .. I . v*  (h  •  COS0  dfi 

jfp  .  v.,.ref.  r'.  x' . r.  .  .  i 

I..  )  cos0^  dSiL 


v,inc  i 


=  i  f  p;<6i.V  =oser  dnr 

•'D 

=  2/"  P^V^r  dMr 

J  0 


(b)  the  hemispherical  reflectivity  (albedo) , 


pv  = 


reflected  flux 


incident  flux 


fWk'1*, in<A)cosei  dt!i 

Q 


_  •'Q 


f  I  .  (ft. )cos0 .  dft . 
/  v,mc  1  11 

•'n 


(c)  the  hemispherical  emissivity, 


e  = 


emitted  flux 


v  black  body  flux  at  T 


=  la 


Jfl  (!2)cos0  dft 

o  V'em 


,BV<Tg> 


-1  A- 

* 


Cu )  JJ  dft 


=  2 


r 


(y)y  dy 


(4. 


(4 
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The  solid-angle  integrals  in  (a)-(c)  are  over  the  upward  hemis¬ 
phere.  Note  that,  by  the  isotropy  of  the  surface,  reallv 

only  depends  on  y^.  Also,  the  above  quantities  are  not  all 
independent.  By  Kirchhoff's  Law  for  surfaces , 

e^(y)  +  P^Cl>)  =  1  (4.18) 

so  that,  if  either  or  is  known,  the  other  is  deter¬ 

mined. 

If  only  the  directional-hemispherical  reflectivity  p^ 
is  given,  nothing  can  be  said  about  the  angular  dependence  of 
the  reflected  radiation.  Therefore  ATRAD  assumes  diffuse  re¬ 
flection,  which  means  p^  in  Ea.  (4.15)  is  independent  of  ur 
and  is  equal  to  p^, 

pJtV’V  =  pv 


If  rather  than  p^  is  known,  Eq.  (4.18)  is  used  to  obtain 

p^,  and  therefore  "p". 

Note  that  the  albedo  p  ,  Eq.  (4.16),  is  not  an  intrinsic 

surface  property  but  depends  on  the  incident  radiation  field 

r  321 

I  ^nc ,  a  fact  familiar  to  experimenters.  J  Only  when  p^  is 
independent  of  angle  does  the  albedo  become  independent  of  in¬ 
tensity,  in  which  case  Eq.  (4.16)  reduces  to 

py  = p;  • 


When  only  an  albedo  is  available,  ATRAD  uses  the  approximation 
that  p"”  is  independent  of  both  yr  and  y^  ,  leading  to 
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To  obtain  the  emissivity,  observe  that  if  the  emission  is 
isotropic  or  if  p^  is  independent  of  angle,  then 


from  Ecr.  (4.18).  Since  the  second  assumption  holds,  e  can 
be  obtained  from  py.  ATRAD  then  invokes  the  assumption  of 
isotropic  emission  to  obtain 


from  Eq.  (4.17). 

A  large  amount  of  albedo  data  is  available , [ 33 ^  particularly 
as  it  relates  to  remote  sensing  from  aircraft  and  satellites.  Much 
more  such  data  remains  classified.  However,  as  yet  no  compre¬ 
hensive  tabulation  of  albedos,  bidirectional  reflectivities,  etc., 
has  been  incorporated  into  ATRAD;  data  for  individual  problems 
will  be  obtained  as  needed. 


4.5  TRANSMISSION  FUNCTIONS 


Line-by-line  absorption  coefficients  have  been  calculated 
and  measured  during  the  last  two  decades  for  all  of  the  important 
atmospheric  absorbers.  However ,  it  is  computationally  infeasible 
to  use  frequency  groups  small  enough  to  resolve  individual  lines 
in  ATRAD.  Therefore,  frequency-averaged  absorption  coefficient 
data  in  the  form  of  transmission  functions  are  employed. 

In  terms  of  the  optical  depth  for  line  absorption 


"v, ab (zi 'z2> 


a  dine) 
v 


(z)  dz , 


(4.19) 


the  transmission  function  is 


70 


SSS-R-72-1255 


TAv  '  2 2  rV) 


_1  f  “Xv,ab(zi'22>/H 
Av  J.  e  dv 

•'Av 


(4.20) 


The  frequency  interval  Av  in  general  includes  many  lines. 

The  computation  of  transmission  functions  has  a  long 
history.  The  earliest  attempts  were  based  on  band  models,  in 
which  simple  analytical  representations  of  line  strengths, 
positions,  and  shapes  were  assumed.  As  accurate  line-by-line 
absorption  data  has  become  available,  both  from  theory  and  ex¬ 
periment,  transmission  function  computations  have  incorporated 
it.  However,  such  detailed  line-by-line  transmission  function 
computations  are  very  expensive  in  terms  of  computer  time.  Con¬ 
sidering  that  integration  steps  as  small  as  0.001  cm"1  may  be 
required,  and  that  the  region  of  significant  absorption  extends 
from  10,000  cm  (ly)  to  250  cm"1  (40y),  with  only  a  few  gaps, 
the  magnitude  of  the  problem  is  apparent.  As  an  example,  Kylet34J 
used  15  minutes  of  GDC  6600  time  to  compute  transmission  func¬ 
tions  between  a  single  pair  of  atmospheric  levels  zl  and  z2, 
for  a  single  value  of  y,  for  the  wavelength  interval  1.7y  -*20y 
with  Av  -  1  cm  .  Multiply  this  by  the  number  of  angles  N  and 
by  the  number  of  pairs  of  levels  ^(N^-l),  and  the  computing 
time  to  obtain  a  complete  set  of  transmission  functions  becomes 
truly  formidable  (for  4  angles  and  20  levels  Kyle  would  require 
190  hours).  While  it  is  foreseeable  that,  with  ingenuity,  such 
computing  times  could  be  reduced  substantially,  it  is  apparent 
that  for  the  present  ATRAD  cannot  afford  to  calculate  transmis¬ 
sion  functions  line-by-line  for  each  new  problem.  Therefore, 
a  less  exact  method  was  sought. 

A  rather  thorough  literature  survey  revealed  numerous 
less  exact  methods  for  computing  transmission  functions.  These 
ranged  from  simple  parameterized  band  models  such  as  Elsasser's 
and  G°°dy's  through  more  sophisticated  band  models  such  as 

Wyatt's,  to  approaches  such  as  those  of  Smith  * 37  ■*  or 
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f  3  8  ] 

McClatchey 1  involving  an  empirically  determined  function  as 
well  as  empirically  determined  parameters.  The  errors  in  the 
methods  of  Smith  or  McClatchey  are  almost  within  the  uncer¬ 
tainties  in  the  line-by-line  calculations  themselves,  and  in¬ 
volve  orders  of  magnitude  less  calculation.  The  method  of 
McClatchey  was  selected  for  inclusion  in  AT PAD  because  of  its 
comprehensive  treatment  of  all  known  atmospheric  absorbers,  its 
accuracy,  and  its  computational  economy  (only  tabular  interpola¬ 
tion  is  required  to  obtain  transmissions) . 

The  method  of  McClatchey,  et.al.,  as  furnished  to  us  in 
subroutine  LOTRAN  (for  LOw-resolution  TRANsmission  functions) , 
contains  separate  values  of  the  transmission  for:  (a)  H20 
vapor;  (b)  the  uniformly  mixed  gases  CC>2 ,  N20,  CH4,  CO,  and  02; 
(c)  03;  (d)  the  H20  8-13y  continuum;  and  (e)  the  N2~N2  4y  con¬ 
tinuum.  In  ATRAD,  the  two  continua  (d)  and  (e)  have  been  placed 
in  subroutines  separate  from  LOTRAN  because  their  transmissions 
are  simple  exponentials  (for  small  enough  Av)  and  they  can  be 
added  to  cont,  the  continuum  part  of  the  absorption  coef¬ 

ficient. 

LOTRAN  furnishes  transmissions  averaged  over  20  cm-1 
intervals  centered  from  350  cm-1  to  50,000  cm-1.  To  obtain 
transmissions  over  larger  wavenumber  intervals,  it  follows 
rigorously  from  the  definition  (4.20)  that  the  20  cm-1  trans¬ 
missions  may  be  averaged;  e.g., 

T400-440  =  2  (T400-420  +  T420-440) 

Hence  ATRAD  may  use  wavenumber  intervals  which  are  any  multiple 
of  20  cm  1,  subject  to  restrictions  on  the  size  of  Av  dis- 
cussed  in  the  prior  semiannual  report.  For  wavenumber  in¬ 
tervals  centered  below  350  cm”1,  the  Goody  random  band  model 
for  the  H20  rotational  band  is  used  in  LOTRAN,  on  the  recom¬ 
mendation  of  one  of  McClatchey' s  colleagues  (J.E.A.  Selby). 
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Parameters  for  this  band  model  were  taken  from  Goody, C 35 ^  Table 

5.5. 


In  order  to  account  for  atmospheric  slant  paths  along 

which  the  pressure  p  and  temperature  T,  and  therefore 
line  ,  . 

•  vary — that  is,  m  order  to  perform  the  integration  in 
Eq.  (4.19) — LOTRAN  employs  a  version  of  the  scaling  (or  effec¬ 
tive  mass)  approximation.  This  approximation  is  an  attempt  to 
sum  up  zi,  z 2  and  the  p-T  variation  of  a^ne(z)  into  a 
single  argument  u  called  the  'effective  absorber  amount.' 
u  is  calculated  using  certain  a  priori  assumptions  concerning 


the  pressure  and  temperature  variation  of  a^lne,  with  parameters 
chosen  to  best  fit  the  real  data.  Following  Goody, t35^  the  scal¬ 
ing  approximation  assumes  the  following  decompo  iition  of 


a 


line 

v 


(2) 


Pa(z)a1  (v)a2  (p,T) 


(the  absorber  density  pa  is  separated  out  because  the  mass 
absorption  coefficient  a1aZ'  an  intrinsic  propertv,  is  preferred 
by  spectroscopists) .  The  stimulated  emission  correction  has  been 
omitted  from  explicit  consideration,  but  the  effects  of  stimu¬ 
lated  emission  are  present  in  the  LOTRAN  model  insofar  as  certain 
parameters  are  determined,  at  least  in  part,  from  measured  trans¬ 
missions*.  The  absorption  optical  depth  (4.19)  then  becomes 

fzz 

Tv,ab(Zl'Z5i)  =  ai  J  Pa  (z)  a2  (p,T)dz 

2  i 

=  kv  u(zirz2)  (4.21) 


*McClatchey,  private  communication. 
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where 

kv  =  ai  (v)a2 (Po /T0) 

"  X/a(Z>  affiff.)  dz  (4-22) 

and  p0  and  T0  are  taken  in  LOTRAN  as  STP  (1  atm. ,  273°K) . 
u  is  the  equivalent  absorber  amount  on  the  vertical  path  be¬ 
tween  zx  and  z2.  The  further  assumption  is  then  made  that 
a2  has  the  form 

a2 (p,T)  =  ApXlTXz 

Different  values  of  yx ,  chosen  to  give  a  best  fit  to  line- 
by-line  transmissions,  are  used  in  LOTRAN  for  each  category 
of  absorber  (a)-(e). 

Unfortunately,  because  of  lack  of  good  data  on  line 
strengths  as  a  function  of  temperature,  it  has  not  been  pos¬ 
sible  to  obtain  best  values  of  y2  i-n  tke  same  way.  An 

attempt  has  been  made  to  include  the  temoerature  dependence 

h 

of  the  line  half-width  by  using  p/x  instead  of  p,  so  that 

at  present  =  “  Yi/2 ,  but  it  is  not  at  all  certain  that 

this  leads  to  less  error  nor  that  T  2  is  indeed  the  correct 

temperature  dependence  of  the  half-width.  A  prescription  for 

putting  temperature  dependence  into  T.  has  been  proposed 

M91  nv 

by  Kondratyev  and  Timofeev. 1  These  authors  claim  that 
transmission  errors  of  up  to  .1  -  .2  may  result  if  tempera¬ 
ture  dependence  is  not  accounted  for,  even  if  one  uses  the 
Curtis-Godson  approximation  instead  of  scaling.  Therefore 
more  work  is  needed  on  this  question. 
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Many  IR  radiative  transfer  studies  employ  the  Curtis- 
Godson  approximation  (see  Armstrong  for  an  excellent  dis¬ 
cussion)  rather  than  the  scaling  approximation.  Very  few 
comparisons  have  been  made,  however,  of  Curtis-Godson  and 
scaling  with  exact  calculations .  Two  such  comparisons  are 
that  of  Kondratyev  and  Timofeev,  discussed  in  the  last  para¬ 
graph,  and  that  of  Zdunkowski  and  Raymond.  The  first  con¬ 

cludes  that  Curtis-Godson  and  scaling  incur  comparable  errors. 
The  second  concludes  that,  while  Curtis-Godson  has  the  edge 
in  accuracy,  the  scaling  approximation  is  accurate  to  at  least 
two  decimal  places  for  a  wide  range  of  path  lengths.  Thus  it 
appears  to  be  justified  to  use  the  scaling  approximation,  es¬ 
pecially  for  applications  such  as  the  present  one  in  which 
errors  in  any  one  frequency  interval  are  mitigated  bv  the  in¬ 
tegration  over  all  frequencies. 


4.6  SEPARATION  OF  INTENSITY  INTO  SOLAR  AND 
DIFFUSE  PARTS 

Because  the  solar  beam  and  the  specular  reflection 
(if  any)  of  the  solar  beam  are  essentially  6-functions  in 
angle,  it  is  customary  before  solving  the  radiative  transfer 
equation  numerically  to  eliminate  these  beams  from  explicit 
consideration  via  the  following  transformation: 


(  I  -  TsPec 
)  v  v 

v  )  y  -  Ysolar 
1v  ~  v 


-1  <  y  <  0 
0  <  y  <  1 


(4.21) 


iv,  the  diffuse  intensity,  may  be  expected  to  behave  'smoothly' 
as  a  function  of  angle,  that  is,  to  be  amenable  to  angular 
quadrature  schemes  which  assume  i  can  be  approximated  by  a 

AV  "]  a  y* 

low-order  polynomial  in  y.  The  solar  intensity  I  and 

the  specularly  reflected  solar  intensity  j^Pec  are  obtained 
by  considering  the  incident  solar  beam  to  be  acted  upon  only 
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by  extinction 


k  and  specular  reflection  p' 
v  v,s 


^solar 

v 


-t  lz)/ji 
e  &  (P  “JJ  0 ) 


jSpsc 

V 


[2tv  (z0)-t  (z)]/y 

pv ,  s  ( I y  I ) e  <S(u+yo) 


where 


Tv(z)  =  f  <v(z')dz' 
•'o 


is  the  optical  depth  measured  from  the  top  of  the  atmosphere. 

If  the  transformation  (4.21)  is  applied  to  the  radia¬ 
tive  transfer  equation  (4.6),  the  boundary  conditions  (4.10) 
and  (4.14),  and  the  flux  equation  (4.7),  the  results  are: 


3i 


9z 


K  X 
V  V 


=  a'B 
v  v 


Pv (z,y,y 1 )iv(z,y* )dy' 


3.  S 


v"v  -  *  _  v  •  ”  r  0 


-t  ,  (z)/y, 


4tt  v 


p  (z,y,y0)  e 


3,,S,>  -[2tv(z0)-tv(z)]  /y0 


+  IT  pi.s(l,o>Vz'M'-u»)e 


(4.22) 


iv(0,y)  =  0, 


0  <  y  <  1 


(4.23) 
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iv(z*'l,)  "  e;<IH)Bv(Tg)  +  p;  (|ji|)iv(z0,| 


f 

*  r\ 


+  2  I  d^’  P"  dllPl,S,)iv(z0,Jl')Jl' 

o  '  '  V  u 


+  IT  >*.  p;,d(ip|-p0)  -6  Tv<‘!,)/|1" 


-1  <  y  <  0 


Vz>  “  2lr  J i  Piv(zfy)dy  + 


(4.24) 


posv  e 


-Tv(z)/y0 


-  y 


°Svpv, 


(y0) 


“[2Tv{zo)-Vz)]/y0 


(4.25) 


These  are  the  equations  which  are  solved  numerically  in  ATRAD. 

There  are  two  functions  in  Eqs.  (4.22)  and  (4.24)  which, 
if  they  were  not  smooth  functions  of  angle,  could  destroy  the 
assumed  smoothness  of  iy  as  a  function  of  y;  they  are  the 
Phase  function  Py  and  the  diffuse  bidirectional  reflectivity 
pv,d‘  The  strongest  peak  in  Py  is  due  to  forward  scattering 
by  aerosol  particles;  therefore,  ATRAD  incorporates  an  approxima¬ 
tion  wherein  this  sharp  peak  is  replaced  by  a  6-function,  leaving  a 
truncated  phase  function  which  is  much  smoother  (see  Section 
4.7.5).  The  same  procedure  can  be  applied  to  a  broadened  re¬ 
flection  peak  in  p"^d,  such  as  might  be  caused  by  a  slightly 
disturbed  ocean  surface.  That  is,  the  peak  will  be  approximated 
by  a  6-f unction  and  lumped  into  the  specular  reflectivity  p* 
leavinq  a  truncated  and  much  smoother  diffuse  bidirectional  re¬ 
flectivity  P"/d(trunc)  for  use  in  the  boundary  condition  (4.24). 
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4.7  SCATTERING  TREATMENT 

4.7.1  Rayleigh  Scattering 

The  Rayleigh  volume  scattering  coefficient, 

o  3  Cnf.-l) 2  /.6+3  p  \ 
fi  =  ®TT _ s  ;  Mn  \  p 

v,R  3  A“*N|  \6“7pl  kT 

and  the  Rayleigh  phase  function  for  unpolarized  light 


■  3U+.p) 

(yJ  =  n 


4  +  2p 


n 


are  taken  from  the  review  article  of  Penndorf. f42]  The  various 
quantities  in  these  equations  are: 


ns  -  index  of  refraction  of  air  at  760  mm  Hg 
and  15°C 

A  =  wavelength 

Ns  ~  numker  density  of  air  molecules  at  760  mm  Hg 
and  15°C  =  2.54743  x  io19  cm"3 

Pn  =  depolarization  factor 

k  =  Boltzmann's  constant 

ys  =  cosine  of  scattering  angle 


ng  is  calculated  according  to  the  new  empirical  formula  of  Peck 
and  Reeder,  which  supersedes  the  older  Edlen  formula.  [42] 
Pn,  which  is  a  measure  of  the  anisotropy  of  the  air  molecules, 
is  taken  as  0.0303  following  Gucker,  et.al.^44^ 
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4.7.2  Mie  Scattering ,  Single  Sphere 

Mie  scattering  of  a  plane  wave  of  wavelength  A  from 
an  aerosol  particle  of  radius  a  is  described  rigorously  by 
the  following  infinite  series: 


°ext  *  57  £  (2n+1>  He(an‘bn1 


n=l 


(4 


sea  2tt 


Y  (2n+X) ( I an| 2  +  |bn|2) 


n=l 


a  ,  =  a  .  -  a 

abs  ext  sea 


(4 


11 


Y  HTHT3T[an(“'n'),Tn(l,s)  +  V^W] 


n=l 


X2 


Y  HTH7lf[an<a'm>Tn(ys)  +  Va'm)Vl,s>] 

n=l 


(4 


(4 


where 


a 


2ira 

A 


m  =  complex  index  of  refraction 


nl  -  in2 


6  =  ma 


.26) 


.27) 


.28) 


.29) 
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^n^01)  ^(3)  "  m  1lJn(3)^(a) 

dn  ”  Cn(a)^(3)  -  m  ^(3)?^  (a) 

^n(3)^(a)  -  m  if»n(a)^(g) 

bn  "  ( 3 )  C'J(a)  -  in  ?n(a)^(0) 

(y)  =  p^  (y) 

Tn (y)  =  yTTn(y)  -  (1-y2)  ir^(y)  . 


The  are  Legendre  polynomials  and  the  ^  are  Ricatti- 

Bessel  functions.  aext  •  asca  •  an(3  °abs  are  the  extinction, 
scattering,  and  absorption  cross  sections  of  the  particle.  i^ 
and  i2  describe  the  patterns  cf  scattered  intensity  polarized 
perpendicular  and  parallel  to  the  scattering  plane  (determined 
by  incident  and  scattered  beams) ,  respectively.  For  the  un¬ 
polarized  light  treatment  assumed  in  ATRAD,  the  pattern  of 
scattered  intensity  is  proportional  to  i1  +  i2. 

It  can  be  shown  that  agca  is  related  to  the  integral 
of  i^  +  i2  over  all  scattering  angles  : 


a 

sea 


(a) 


[il(a,ys)  +  i2 (a,ys) ]d^s 


(4.30) 


The  numerical  evaluation  of  the  infinite  series,  Eqs . 
(4.26)  through  (4.29),  while  seemingly  straightforward,  has 
sources  of  error  which  are  related  to  the  accumulation  of 
round-off  error  in  recursively  calculating  ,  5  ,  tt  , 

an(3  .  Techniques  to  eliminate  these  errors  have  been 
discussed  by  numerous  authors.  Dave,^45^  Kattawar  and 
Plass ,  f  46  ]  ancjj  Denman,  Heller,  and  Pangonis^4^  are  the  pri¬ 
mary  sources  of  the  present  ATRAD  formulation  of  the  series 
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summation,  in  subroutine  MIESCT.  Calculated  results  were  com¬ 
pared  with  che  graphs  of  Kerker^15^  and  the  tables  of  Deirmend' 
jian^^  and  of  Denman,  et.al. ^ ^  Agreement  was  found  in  all 
cases  with  the  published  results. 

The  Mie  scattering  calculation  requires  the  aerosol 
complex  index  of  refraction  mv(z)  as  a  function  of  frequency 
and  height.  This  is  the  least  well-known,  for  the  real  atmos¬ 
phere,  of  all  the  required  aerosol  data.  Excellent  m^  data 
are  available  for  water;  we  have  selected  the  measurements  of 
Irvine  and  Pollack, Rusk,  et.al. , and  Robertson  and 
Williams for  inclusion  in  ATRAD  (subroutine  IOR) .  Even 
for  water,  however,  the  uncertainty  in  the  real  part  of  m^ 
is  as  high  as  5  percent  and  in  the  imaginary  part  as  high  as 

20- 25  percent.  Somewhat  less  accurate  data  are  available 

•  r  49 ]  (521 

for  ice,  aqueous  solutions  of  NaCl,  1  carbonaceous  aero- 

[53  54]  [55i 

sol,  '  and  quartz,  all  of  which  have  been  included 

in  ATRAD.  The  above  substances  have  each  been  studied  by 
several  investigators,  so  at  least  the  data  have  been  measured 
independently.  No  such  independent  evaluation  exists  for  the 
index  of  refraction  measurements  of  Volz,  ^6/57]  who  studied 
sea  salt,  dust,  Sahara  dust,  and  several  samples  of  water- 
soluble  aerosol  material;  nevertheless,  these  data  have  also 
been  included  in  ATRAD. 

Aerosol  materials,  of  course,  have  highly  variable 
compositions  and  are  subject  to  uncertainties  regarding  their 
chemical  and  physical  states. 

A  large  number  of  transmission  spectra  have  been  taken 
for  atmospheric  dust  and  for  individual  components  of  atmos¬ 
pheric  dust.  ,60]  These  c;pectra  indicate  strong  dust 

absorption  in  the  atmospheric  window  region  8-13y,  which  could 
have  a  pronounced  effect  on  radiative  heating.  However, 
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more  data  than  the  transmission  spectrum  of  a  substance  are 
needed  to  deduce  the  index  of  refraction.  Even  less  is  known 
about  the  index  of  refraction  of  many  other  inorganic  aerosol 
constituents.  As  for  the  organic  aerosol  fraction,  which  is 
considerable,  ^  the  lack  of  data  is  almost  complete. 

4 • 7 . 3  Mie  Scattering,  Polydispersion 

In  reality,  the  atmospheric  aerosol  consists  of 
particles  with  a  wide  range  of  radii  a.  This  variation  is 
described  by  a  size  distribution  function. 


Naer(z)  n(a,z) 


a  .  <  a  <  a 

mm  —  —  max 


where  Naer(z)  is  the  total  number  density  of  aerosol  particles 
at  height  z  and  n(a,z)  is  the  normalized  distribution  of 
radii  at  height  z,  defined  such  that 


n (a , z) da  =  fraction  of  particles  at  height  z 
with  radii  in  (a,a+da) 


f 

•'a 


max 


n (a, z) da 


min 


The  separation  of  ^aer  from  n  is  convenient  because  in 
practical  applications  we  eften  assume  n  is  independent  of 
z  .  The  volume  scattering  and  absorption  coefficients  for 
an  aerosol  described  by  such  a  size  distribution  are 


3v,M(z)  “  Naer(z)  f  ^  asca(a)  n(a'z)da 

mm 


(4.31) 
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Vm(z)  =  "aer'2*  /  ”'aX  ®abs(ci)  n(a,Z)da  .  (4.32) 

mm 

The  pattern  of  scattered  intensity  will  be  proportional  to 

r  a 

Naer^  J  maX  +  i2^'a,ys^  n(a,z)da  . 

a.  • 
mm 


Normalizing  this  so  that  its  integral  over  all  scattering  direc¬ 
tions  is  4tt  ,  we  obtain  the  Mie  phase  function 


/a 

maX  n(a,z)da 


pv,j.‘z-«.>  =  4" 


a  . 
mm 


f  f  max 

J ^  dns  J  fi].(a'  ys)  +  i2(a,ys))  n  (a,  z)  da 


mm 


Inverting  the  order  of  integration  in  the  denominator  and 
employing  the  relationships  of  Eqs.  (4.30)  and  (4.31), 


A2N  (z) 

P  (z  u  )  =  aer 
V,M<  'V  ^VfM(z) 


fa 

J  max  [i1(a,ys)  +  i2(a,ys)]  n(a,z)da 


mm 


(4.33) 


Various  analytic  forms  for  n(a,z)  have  been  proposed. 
Perhaps  the  most  popular  is  the  Junge  power-law  distribution^62^ 
(ca"’);  in  the  American  literature  Derimendjian ' s  modified 
gamma  distribution [63]  (caY  e-ba6)  for  haze  and  cloud  is  also 
used  extensively.  Log-normal  distributions 
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/  1  -In2  (r/rQ)/2o2\ 

a  /2tt  a  J 

are  used  by  some  European  and  Russian  aerosol  researchers . [64 ' 65] 
Twomey ,  et.al.,[66]  employ  Gaussian  distributions,  but  these  are 
not  in  general  use  because  most  available  data  indicate  an 
asymmetrical  distribution  of  aerosol  sizes.  With  the  exception 
of  the  Junge  distribution,  which  is  unrealistically  monotonic 
for  all  a,  the  data  do  not  favor  one  size  distribution  formula 
over  another.  All  of  the  above  are  included  as  options  in 
ATRAD. 


The  dependence  of  Naer(z)  on  z  for  the  cloudless 
atmosphere  has  been  investigated  extensively ;  ^7'  ^8  •  69  >  70 , 71] 
certain  regularities  have  been  determined.  Eltermann,  *-72-*  for 
example,  has  proposed  a  clear  standard  atmosphere  with  a  pre¬ 
scription  for  Naer(z)  which  has  been  included  in  ATRAD  as 
an  option.  In  clouds,  Naer  may  vary  widely, [64]  but  most 
models  ignore  this  and  assume  homogeneous  clouds.  The  varia¬ 
tion  of  Naer  (as  well  as  n)  near  a  cloud's  upper  and  lower 
boundaries  may  have  important  consequences  for  atmospheric 
heating,  however,  and  should  not  be  ignored.  Such  variations 
can  be  incorporated  quite  straightforwardly,  within  the  frame¬ 
work  of  the  ATRAD  formalism. 

A  further  complication  for  reu  stic  aerosols  is  that 
the  size  distribution  and  the  index  of  refraction  depend  on 
the  relative  humidity . ^73 • 74 ^  Hygroscopic  aerosols  may  begin 
to  take  up  water  at  a  relative  humidity  as  small  as  50  percent. 
Hence  it  is  impossible  to  specify  realistic  atmospheric  struc¬ 
ture  data  for  ATRAD  without  taking  this  variable  into  considera 
tion. 
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4 • 7 • 4  Numerical  Integration  Over  Aerosol  Size  Distribution 

Techniques  for  performing  the  numerical  integrations 
over  size  distribution  in  Eqs.  (4.31)  —  (4.33)  are  far  from 
standard.  Most  investigators  use  trapezoidal  quadrature  on 
the  grounds  that  osca  ,  0abs  ,  ^  ,  and  i2  oscillate  so 

rapidly  and  non-uniformly  with  a  that  more  sophisticated 
quadrature  schemes  are  less  accurate.  •  This,  of  course,  is  a 
tacit  admission  that  their  integration  meshes  are  not  fine 
enough  to  resolve  the  oscillations.  However,  the  present 
authors  performed  several  integrations  over  size  in  which 
the  oscillations  we re  resolved,  and  the  improvement  in  accu¬ 
racy  was  too  small  (never  more  than  a  few  percent)  to  warrant 
the  additional  computing  cost  —  especially  in  light  of  uncer¬ 
tainties  in  the  Mie  input  data  for  a  real  aerosol.  Therefore, 
Romberg  quadrature,  which  is  based  on  multiple  trapezoidal 
quadratures,  was  selected  for  use  in  ATRAD. 

The  comprehensive  study  of  Deirmendjian  ^63  ■*  proved 
useful  in  ascertaining  suita' ..'.t1  integration  increments  Act 
for  the  Romberg  quadrature.  However,  Deirmendjian '  s  Act's 
were  chosen  by  trial-and-error ;  he  did  not  determine  a  univer¬ 
sal  prescription  for  selecting  Act's  small  enough  to  ensure 

accuracy  yet  large  enough  to  avoid  wasting  computer  time. 

f 76  771  1 

Studies  by  Dave  '  indicate  the  danger  of  selecting  Act 

too  large  and  show  that  a  Act  which  is  satisfactory  for  inte¬ 
grating  and  m  ,  and  M  for  small  angles,  may 

be  unsatisfactory  for  integrating  ^  for  larger  angles. 
However,  Dave,  whose  work  was  restricted  to  water  droplets  in 
the  visible  spectrum,  provides  no  rule  for  automatic  selection 
of  Act  . 

We  have  devised  an  automatic  selection  scheme  for  Act 
which  has  proven  successful  in  duplicating  the  extensive 
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T6  3  1 

tables  of  Deirmendjian . 1  J  It  is  based  on  our  observations 
that : 


(a)  as  the  upper  limit  of  integration  a 
increases,  computing  expense  increases 
inordinately,  because  the  number  of 

terms  needed  in  the  Mie  series  [Egs.  (4.26) 
through  (4.29)]  is  proportional  to  a  ; 

(b)  the  number  of  points  necessary  in  the 

integration  mesh  increases  more  or  less 

linearly  with  a  •  ; 

1  max 

(c)  the  angular  interval  A6  needed  to  re¬ 
solve  the  forward  peak  in  P  is  in- 

v  ,m 

versely  proportional  to  amax  •  but  away 
from  the  peak  progressively  larger  angular 
steps  may  be  taken . 

In  light  of  observations  (a)  and  (b) ,  it  was  decided 
to  perform  the  crudest  of  the  Romberg  trapezoidal  integrations 
with  a  step  size 


Aa  = 
o 


a  -  a  . 
max  mm 


N. 


a 


where 
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N  = 
a 


(Nmir 

/  M 


\N  .  +  B(a  -N  .  ) 

\  min  max  min7 

l 

\  max 


a 

max 

N  . 
min 

a 

max 


<  N  . 

—  min 


<  a 

max 


>  N  . 
min 


N  -N  . 

<  M  +  max  min 
—  min  B 

Nmax  Nmin 
+  - - - 


and  to  terminate  this  integration  (and‘,  therefore,  the  more 

refined  integrations  with  Aa  =  ^Act  ,  Act  =  ^Aa  ,  etc.) 

AO  4  O 

whenever  the  relative  changes  in  3  ,,  ,  a  ,  and  P 

^  V  ,M  '  v,M  '  V,M 

resulting  from  addition  of  the  next  term  are  less  than  6 

o 

The  current  ATRAD  values  of  these  parameters  are 


N  . 
min 


25 


N 

max 


80 


B  =  0.55  6 

o 


0.001 


which  are  sufficient  to  reproduce  the  tables  of  Deirmendjian. 

By  terminating  the  integration  in  this  fashion,  unnecessary 
computation  is  avoided. 

The  forward  (diffraction)  peak  in  P  for  a  particle 

_  .  I  1 v 

of  size  parameter  a  is  proportional  to1  1 


F(z) 


2  J^zjj 


2 


where  z  -  asinG  ,  and  is  the  Bessel  function  of  order 

one.  This  quantity  is  shown  in  Figure  4.3. 

The  actual  forward  peak  will  be  a  composite  of  the 
peaks  for  all  ae t amin' amax^ '  T^e  s^arPest  of  the  component 
peaks  is  due  to  a  =  a  ,  and  therefore  in  order  to  resolve 

lltciX 

this  peak,  which  varies  substantially  over  Az  =  a  AG  =  1  . 

max  ’ 

we  pick  an  angular  increment 
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Figure  4.3  —  Graph  of  F(z) 


A6  =  ,5 - i 

o  2  a 


Actually,  because  the  t^'s  and  xn ' s  entering  the  Mie 
series  [Eqs.  (4.28)  and  (4.29)]  are  pre-calculated  with  a 
predetermined  angular  mesh  A0fcab  ,  A0q  is  selected  to  be 
the  closest  multiple  of  A0tab  to  1/2  ctmax  .  A0q  is  also 
limited  to  be  less  than  0^  (currently  1°).  The  angular  mesh 
on  which  P^M  is  calculated  is  then  taken  as  a  user-supplied 
(see  NSTEP ,  Appendix  A)  and  even  number  of  steps  of  A0  , 
then  2A0q  ,  then  4A0q  ,  etc. ,  with  the  restriction  that  no 
step  is  larger  than  0max  (currently  2.5°).  The  number  of 
steps  of  each  size  is  even  so  that  a  piecewise  Simpson  quad¬ 
rature  can  be  used  to  integrate  P  in  case  truncation 
.  v,M 

(Section  4.7.5)  or  renormalization  (Section  4.7.7)  is 

necessary. 
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4.7.5  Phase  Function  Truncation 

Hansen^  ^and  Potter  have  shown  that  the  forward 

peak  of  the  phase  function  Py,M  can  be  truncated  without 

significantly  changing  the  results  of  the  transport  calculation. 

Thus,  for  6 e  £ 0  r  0  ],  where  0t  is  the  truncation  angle,  the 

real  phase  function  Pv(z,0)  is  'replaced'  (to  within  a 

normalization  factor)  by  a  truncated  phase  function  P*"  (z  0) 

.  v  m  f  ' 

m  the  manner  of  Figure  4.4.  The  logarithm  of  P^  is  taken 

to  be  a  linear  function  of  0  on  [0 , ©fc]  and  the' values  of 

Pv,M  its  derivative  are  matched  to  those  of  P  at  0  =  0^ 

'  v,M  * 

Neither  Hansen  nor  Potter  propose  a  general  method  for 

selecting  0  .  Therefore  we  devised  the  following  procedure: 

If  Pv,M  for  0=0  is  less  than  a  prescribed  P  (currently 

'  '  v,M  1s  truncated.  Otherwise,  truncation  is  attempted 

at  a  succession  of  O^s  beginning  near  0  =  0  until  either  a 

successful  truncation  is  achieved  (such  that  Pfc  (z,0)  <  P  ) 

or  0fc  exceeds  0T  (currently  20° )  .  If  0fc  exceeds  0T  ""first, 

Pv,M  ls  truncated  at  0  without  regard  to  P  ,  provided 

only  that  pv,M(Z/0)  <  pV/M^z/°).  It  the  latter  procedure  is 

unsuccessful,  P^  M  is  not  truncated. 

Mathematically,  trunction  amounts  to  approximating  the 
forward  peak,  P^M  -  P^M,  by  a  6-function: 

PV,M  “  (PV,M  "  PV,M)  +  PV,M 

-  477F  6  (h '  -  ?2)  +  pfc 

V,M 


where  F  may  be  determined  by  integrating  both  sides  over  all 
scattering  angles 
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F  =  W  jT  (Pv,M  - 

'  1  /*Tr  4- 

“2  *Pv,M^2'0*  ~  pV/MCz/0)Jsin6  d6 

1  /*et 

7  Jq  fPv,M(z'0)  "  Pv/M(2/6^sin6  de  (4.34) 


F  is  determined  numerically  (in  subroutine  TRUNCT)  using  a 
Simpson  quadrature. 

Substituting  the  6-function  approximation  into  the 
scattering  terms  of  the  radiative  transfer  equation  (4.6), 

and  eliminating  subscripts  and  unnecessary  arguments  for 
brevity. 


scat,  terms  =  p  (fi/fi*  )I  ($•  )dft*  -  0  I  (ti) 

“  I?  y*[4TTF  6(fl»-^)  +  pt]i(^')dr;:'  -  0  I(tj) 

=  3  (i-F )  r  pfc  $  , ti ' )  _  . . .  _ . 

4  7T  J  1-F  P  (^  )dfi '  -  0(l-F)I(f2) 

=  y*p'  (^,t2')Kfi»)d«'  -  e'Kfi)  , 


where 
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Therefore,  the  scattering  terms  have  their  original  form,  but 
P  and  3  are  replaced  by  P'  and.  3'  defined  above. 


4.7.6  Azimuthal  Integration 

The  aximuthally-averaged  phase  function  P  (Eq.  (4.5)), 
not  Pv  itself,  is  required  in  the  ATRAD  formulation  of  radia¬ 
tive  transfer.  The  azimuthal  average  of  the  Rayleigh  phase 
function  (Section  4.7.1)  is  computed  analytically  in  ATRAD,  from 


C,  -2  TT 

PR(y,u')  =  —  '  "  -  -  "2 


1  f 

2ir  / 
•'0 


(1  +  c,  y2)d(J> 


s 


2^  f  tl+c2  (uu1 )  2+2c2py'  /  (1-y  2 )  (l~y  '  2 )  COs<|> 

Jo 

+  c2  (1-M2)  (1-p  1  2  )  cos2(j)]  d(J) 


=  Cj  [1  +  c2 (yy ' ) 2  +  j  c2 (1-y 2 ,  (1-y'2)] 


where 


3(1  +  pn)  _  1  -  pn 

C1  -  "4  +  2pn  °2  =  1  +  Pn 

The  aziiituthal  average  of  the  Mie  phase  function  is 
computed  numerically.  The  techniques  are  described  in  the 
prior  semi-annual  report  of  this  contract. 1  J  The  only  change 
is  that  the  quadrature  has  been  changed  from  trapezoidal  to 
Simpson  in  order  to  reduce  errors  affecting  the  renormalization 
procedure  (section  4.7.7). 
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4.7.7  Renormalization  of  the  Phase  Function 

There  are  two  types  of  renormalization  of  the  Mie  phase 

function  in  ATRAD.  The  first  operates  on  the  Mie  phase  function 

P^  M  before  azimuthal  averaging,  and  consists  in  multiplying 

all  values  of  P  by  the  renormalization  factor 

v  ,M  J 


2 


P 


V  ,M 


t*,Us)dyB 


(4.35) 


where  the  integral  in  the  denominator  is  computed  numerically 
using  Simpson's  Rule. 

The  second  type  of  renormalization  operates  on  the 

aximuthally  averaged  phase  function  P  M  and  is  required  by 

T  2  1 

the  Grant  and  Hunt  algorithm  in  order  to  preserve  flux 
conservation.  In  terms  of  the  Gaussian  angular  quadrature 
mesh  0  <  y x  <  ...  <  ym  £  1  defined  in  the  previous  semiannual,  ^ 
the  numerical  equivalent  of  the  normalization  condition 


(z,y,y' )dy 


2 


(4.36) 


is 

m 

y]  Cj  Pjy.  =  2  (k  =  l,...,m)  (4.37) 

j=l 

where 
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Cz,jij  ,yk) 


There  are  several  reasons  why  P^  M  will  in  general  fail  to 
satisfy  Eq.  (4.37).  One  is  that  Gaussian  quadrature  does  not 
accurately  represent  Eq.  (4.36);  that  is,  M  may  not  closely 
approximate  a  polynomial  in  y.  Also,  errors  are  introduced  in 
each  step  of  the  process  by  which  Pv  M  is  constructed: 

(a)  the  numerical  computation  of  the  Mie 
series  (4.26)- (4.29) 

(b)  the  numerical  integration  over  the 
aerosol  size  distribution,  Eq.  (4.33) 

(c)  the  numerical  integration  to  obtain 
the  truncation  factor  F,  Eq.  (4.34) 

(d)  the  numerical  integration  over  azimuth. 


Tests  have  been  run  to  determine  ATRAD's  errors  from  each  of 
the  four  sources  (a)-(d);  it  has  been  determined  that  in  no 
case  do  they  exceed  i%.  Hence,  the  necessity  for  renormaliza¬ 
tion  will  be  almost  entirely  due  to  the  inadequacy  of  Eq.  (4.37) 
as  an  approximation  to  Eq.  (4.36). 


ing  Eq 
calc 


Let  us  assume  that  the  corrected  values  satisfy 

(4.37),  are  to  be  obtained  from  the  calculated  values 
by  applying  a  multiplicative  correction. 


P 


jk 


(1  +  e 


calc 

pjk 


94 


41 


SSS-R-72-1255 


Unfortunately,  there  are  ■■■■  w“  Jm>  unknowns  e. 


— —  —  unknowns  (the  number  is 

reduced  from  m2  by  the  symmetry  requirement  e,  .  =  e.,  )  and 
there  arf  only  m  equations  (Eqs.  (4.37))  to  determine  them. 
This  underspecif ication  can  be  handled  in  many  ways,  no  one  of 
which  seems  much  preferable  on  physical  grounds  over  ary  other. 
Grant*  corrects  only  the  diagonal  elements,  6^,  and 

so  has  a  determinate  set  of  equations  for  the  e . .  He  points 


out  that  the  phase  matrices  p 


jk 


are  usually  strongly  diagonally 


dominant  and  that  his  procedure  thus  incorporates  the  entire 
correction  where  it  makes  the  smallest  relative  change. 


An  alternative  procedure  has  been  investigated.  Assume 


'  jk 


ej  +  ek- 


Then  Eqs.  (4.37)  become 

m 

E  (cjpjk  +  Dk6jk)ej  =  2  "  bk  (k-l,...,m) 
j=l 


where 


calc 

Pjk 


and  is  the  Kronecker  delta.  This  is  a  determinate  set  of 

linear  equations  for  the  Ej. 

Both  Grant's  method  and  our  method  of  renormalization 
have  been  included  as  options  in  ATRAD.  The  effect  of  our 


♦Private  communication. 
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method  is,  in  practice,  similar  to  Grant's;  the  largest  cor¬ 
rections  are  made  to  the  diagonal  elements  of  Pj^*  Progressively 
smaller  relative  corrections  result  when  an  element  is  farther 
from  the  diagonal.  However,  our  method  offers  the  advantage  of 
spreading  the  correction,  while  Grant's  procedure  concentrates 
it  entirely  on  the  diagonal  element.  In  our  opinion  it  is 
difficult  to  justify  correcting  only  a  subset  of  the  p^'s. 

Several  test  calculations  have  been  run  in  which  the 
Henyey* Greenstein  phase  function, 

1  -  g2 _ 

(1  +  g2  -  2gyrf> 3/2 


was  used  for  P  ,  with  g  =  0.8  and  g  =  0.95  so  that 
large  forward  scattering  peaks  were  being  tested.  Twelve 
Gaussian  angles  (m=6)  were  used.  For  these  examples,  it  was 
found  that  reasonably  large  (up  to  8%  for  g=0.8  and  up  to  20% 
for  g=0.95)  renormalization  corrections  were  necessary.  The 
results  point  up  to  the  importance  of  truncation  to  ensure 
that  P  „  can  be  treated  by  low-order  Gaussian  quadrature. 


4.7.8  Total  Phase  Function 


The  total  phase  function  P  is  obtained  from  P 


R 


and  P^  M  according  to 


^  ev,R  PR  +  0V,M  Pv,M 

V  8.. 


where 


^V,R  +  0V,M 


is  the  total  volume  scattering  coefficient. 
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4.8  NUMERICAL  SOLUTION  OF  THE  RADIATIVE  TRANSFER 
PROBLEM 

The  numerical  methods  employed  to  solve  the  radiative 

transfer  problem,  Eqs.  (.4 . 22)  -  (4 . 25) ,  were  presented  in  the 

[21 

previous  semiannual.  Therefore  only  changes,  improvements, 
and  items  omitted  from  that  report  are  discussed  in  this  section 

4.8.1  Source  Doubling 

The  Planck  source  doubling  scheme  has  been  derived  in 
a  computationally  more  efficient  form.  The  details  are  given 
in  the  paper  we  presented  at  the  AMS  Conference  on  Atmospheric 
Radiation  (Fort  Collins,  Colo.,  Aug  1972);  that  paper  is  re¬ 
produced  in  Appendix  C. 


4.8.2  Treatment  of  Overlapping  Bands 

In  many  frequency  groups  Av,  more  than  one  atmospheric 
constituent  may  contribute  to  the  absorption.  ATRAD  then  fits 
the  transmission  function  T^"^  of  each  individual  constituent 
with  an  exponential  sum. 


Tinv)  <»> 


M  ' 

z 

•1=1 


-)Jn)u(n) 


Superscript  1  refers  to  H2O,  superscript  2  to  the  uniformly 
mixed  gases  CC^,  etc.,  and  superscript  3  to  0^.  The  product  of 
these  individual  transmission  functions  is  taken  as  the  total 
transmission. 


i=i 


z 


M 


(3) 


=1 


a  (1)  _  (2)  (3  ) 

ai  aj  a* 
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Then  M  =  monochromatic  problems  are  calculated, 
one  for  each  term  of  the  triple  sum,  in  each,  of  which  the  opti 
cal  depth  of  zone  Iz,z']  is 


At(z,z')  =  kf^u^  (z,z')  +  kj^u  ^  (z,z') 

+  k,(3)u(3)  (  Z  ,  Z  1  ) 


and  the  line  absorption  coefficient  at  z  is 

line ,  .  3*/  .  \ 

a  (z)  =  -  ^  At(z,z') 


k.'1)  3 

l 


3z 


u(1)  (z,z ' ) 


kj 


(2) 


3_ 

3z 


U(2)  (z,z' ) 


-  k 


(3)  8_ 
i  dz 


u(3)  (z,z '  ) 


The  z-derivatives  of  the  u's  are  simply  the  integrands  in  the 
expressions  (4.22)  defining  the  u's. 

4.8.3  Special  Case  of  No  Scattering  in  the  Grant-Hunt 
Algorithm 

In  the  case  of  negligible  scattering  it  is  not  neces¬ 
sary  to  build  up  the  solution  from  optically  thin  layers; 
simple  analytic  expressions  for  the  reflection  and  transmis¬ 
sion  matrices  and  Planck  source  vectors  can  be  derived  for 
zones  of  arbitrary  optical  depth.  These  analytic  expressions 
are  obtained  from  the  Grant-Hunt  algorithm  by  taking  the  limit 
as  the  primary-layer  thickness  Aip  goes  to  zero.  For  the 
reflection  and  transmission  matrices  of  a  zone  of  optical 
thickness  Ax  ,  we  obtain 
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For  a  Planck  source  which  varies  linearly  in  x, 


VT>  =  VTo>  +  Bi(T"T*) 


where 


B'  = 
v 


Bv(T,tAT).  -  Bv  (T,,) 

At 


we  obtain  for  the  component  of  each  source  vector. 


ZPlkj  =Bv(To)(1_e 


-Ax/y . 


)  + 


At  -Ax/y.  ] 

—  -  (1  -  e  3) 

yj 


EPlk.  "  Bv(t.)(1  -  e 


-Ax/yj 


)  +  ^ 


i  -  a  +  e'iT/u’ 


These  special  case  results  are  incorporated  in  ATRAD. 


4.8.4  Exponential- Sum  F i tt ing 

The  method  of  Cantor  and  Evans ^  ^  for  exponential- 

sum  fitting  is  now  incorporated  into  ATRAD,  although  numerous 
modifications  and  simplifications  have  been  introduced  in 
adapting  the  method  to  the  fitting  of  transmission  functions. 
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The  method  is  briefly  outlined  in  the  following  paragraphs. 

Mb  desire  to  fit  a  transmission  function  T^^(u), 
whose  values  TAv(un)  are  given  at  a  set  of  points 

0  =  Uj  <Uj,  <  ...  ^  Ujj  f 

with  a  positively-weighted  sum  of  exponentials, 

M  -k.u 

TAv<u)  “  E  aie  '  ai  ’  0 

i=l 


The  fit  is  to  be  best  in  the  least-squares  sense,  that  is, 
it  is  to  be  such  that  the  residual 

N 

«  -  2>n 

n=l 

is  minimized.  The  quantities  WR  are  positive  weights  which 
may  be  assigned  arbitrarily.  In  order  to  ensure  that  small 
and  large  values  of  the  transmission  receive  equal  treatment, 
the  weights  have  been  taken  as 

Wn  = 

so  that  R  is  in  fact  a  sum  of  squared  relative  errors. 

The  u  1  s  are  selected  to  be  equally  spaced, 


tAv<V 


-E 


a .  e 

l 


"KiUn 


u  =  (n-l)Au 
n 
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(this  requirement  will  be  partially  relaxed  in  the  near  future) . 
Then  R  may  be  written 


R 


tun> 


H 


E 

i=l 


-  fln-l 
aiei 


where 


-k .  Au 
9i  =  e  1 

Let  us  presume  the  following  ordering  of  the  0^'s: 

O<0.<6  <  ...  <  0  <  1 

—  1  2  M  — 

The  problem  of  minimizing  R  is  partly  linear  (finding 
a i ' • • • ' aM )  and  Partly  nonlinear  (finding  0lf...,0M).  The  algo¬ 
rithm  which  we  shall  describe  below  treats  the  linear  and  non¬ 
linear  parts  of  the  problem  separately,  although  we  shall  see 
that  there  is  some  coupling. 

Given  a  set  of  we  shall  call  the  procedure 

of  finding  the  a1,...,a^  which  minimize  R  with  no  constraints 
"solving  the  linear  problem."  Since  the  a^s  are  unconstrained 
solving  the  linear  problem  may  produce  one  or  more  negative  a^'s 

The  extreme  (and  well-documented)  ill-conditioning  of 
the  linear  problem  is  side-stepped,  as  far  as  the  calculation 
of  the  residual  is  concerned,  by  the  use  of  divided  differences, 
which  is  made  possible  by  the  use  of  equally-spaced  un's  . 

During  the  iteration  between  linear  and  non-linear  problems, 
the  a^  are  of  importance  only  in  determining  which  exponential 
factors  0i  are  to  be  dropped  because  of  a  negative  coefficient. 
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To  find  the  best  0^'s  we  proceed  iteratively.  Imagine 
that  to  the  approximant 


£ 


a .  eT1 

r  i 


0 


we  add  a  new  term  with  exponential  factor  8  and  coefficient 
a.  We  want  to  see  how  the  residual  changes  as  we  increase  s.a 
from  zero;  therefore,  we  evaluate  the  slope  of  R  at  a=0: 


N 

m 

2} 

d 

da 

£  w„kv<v 

a-ef1  -  ae"-1 

n=l  1 

i=l 

) 

JL  r 

m 

a=0 


-  -  2  £  Kn  Tiv(V  -  £  aier1 


n=l 


i=l 


,n~l 


(4.38) 


which  is  simply  a  polynomial  in  0.  Since  we  want  to  minimize 

the  sum  of  squares,  we  pick  0  so  that  this  expression  is  as 

negative  as  possible  (if  it  is  never  negative  for  0  £  0  <_  1, 

we  are  already  at  a  best  fit).  In  other  words,  we  find  the 

absolute  minimum  of  the  polynomial  in  (4.38)  over  0  <_  0  £  1. 

This  search  was  originally  done  by  simply  partitioning  [0,1] 

into  a  large  number  of  points  (^  1,000)  and  finding  the  absolute 

minimum  over  that  set  of  points.  This  kind  of  search  proved 

to  be  quite  expensive  computationally,  however,  since  it 

s  t 

required  evaluating  an  (N-l)  degree  polynomial  (N  ^ 

15-50)  at  more  than  1,000  points,  for  each  iteration  of  the 
fitting  process,  for  each  molecular  species,  in  every  fre¬ 
quency  level.  Therefore,  the  search  has  been  improved 
by  taking  into  account  the  third  derivative  of  the  polynomial. 
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With  this  scheme  large  fractions  of  [0,1]  can  be  eliminated 
from  consideration  and  convergence  is  accelerated  in  the 
remainder  of  the  interval .  The  new  procedure  requires  far 
fewer  polynomial  evaluations. 

When  the  new  exponential  factor  0  is  added  to  the 
previous  ones  and  the  linear  problem  for  the  coefficients  is 
solved  we  obtain  an  approximant 


m 


i=l 


Because  the  polynomial  (4.38)  was  negative  at  0,  we  know 
that  c  >  0;  but  it  may  happen  that  some  of  the  b^'s  are 
non-positive.  In  that  case,  we  move  on  a  straight  line  in 
coefficient  space  from  (a^ , . . . ,a  ,0) f  the  point  represent¬ 
ing  the  old  coefficients,  to  (bt , . . .  ,b  ,c)  ,  the  point  repre¬ 
senting  the  new  coefficients.  We  stop  when  we  encounter  the 
first  zero  coefficient  and  eliminate  the  corresponding  exponential 
factor.  Then  we  re-solve  the  linear  problem  and  repeat  the 
above  process  if  any  of  the  resulting  coefficients  are  non¬ 
positive.  When  solving  the  linear  problem  finally  yields  all 
positive  coefficients,  we  are  ready  to  solve  the  nonlinear 
problem  again. 

In  moving  toward  a  solution  of  the  linear  problem,  we 
constantly  diminish  the  residual  as  we  proceed. 

A  diagram  summarizing  the  algorithm  follows: 
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Start  with  approximant  =  0 

>  Add  new  exponential  factor  6e[0fl] 
which  minimizes  polynomial  (4.38) 

1 

Solve  linear  problem  for  best  coef¬ 
ficients  in  least  squares  sense, 
eliminating  exponential  factors  as 

needed  to  obtain  positive  coefficients. 
- 1 


The  algorithm  is  terminated  when  the  residual  R  increases 
from  the  preceding  iteration  or  when  its  relative  decrease  is 
at  the  noise  level.  Either  situation  indicates  that  a  best 
fit  has  been  obtained  and  that  the  algorithm  is  "spinning  its 
wheels"  due  to  round-off. 

For  the  purposes  of  the  radiation  code,  two  clean-up 
steps  have  been  added  to  the  algorithm  for  the  sake  of  com¬ 
putational  efficiency. 

First,  it  is  often  the  case  that  the  exponential  factors 
0^  occur  in  very  close  pairs.  Since  every  extra  0.  means 
extra  monochromatic  problems  to  solve,  it  seemed  reasonable  to 
replace  each  close  pair  (0i/0i+1)  bY  a  single  with  a 

new  coefficient.  This  procedure  causes  essentially  no  loss 
in  accuracy  in  the  radiation  computation. 

Secondly,  all  coefficients  which  are  small,  a^  £  e, 
have  the  corresponding  terms  dropped  from  the  approximant. 

The  remaining  coefficients  are  renormalized  so  that  they  sum 
to  unity,  in  order  to  preserve  flux  conservation.  Currently 
we  take  e  =  .005,  in  view  of  inaccuracies  in  our  transmission 
functions  of  roughly  1%. 

The  polynomial  (4.38)  plays  an  important  role  in  this 
theory.  Based  on  its  minimum  value  for  a  given  approximant. 
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a  bound  can  be  derived  on  the  "distance"  (in  the  least  squares 
sense)  between  the  best  approximant  and  the  given  approximant. 
Also,  a  lower  bound  on  the  amount  of  improvement  that  will  re¬ 
sult  from  the  next  iterative  step  can  be  derived.  Taken  together, 
these  bounds  guarantee  the  convergence  of  the  algorithm. 

If  the  best  fit  is  not  perfect,  then  it  can  be  shown 
that  the  exponential  factor.*  and  coefficients  of  the  best 
approximant  are  unique.  In  fact,  the  exponential  factors  will 
all  be  from  among  the  roots  between  zero  and  one  of  the  final 
polynomial  (4.38). 

4.9  ATRAD  FLOW  DIAGRAM 

A  simplified  flow  diagram  of  ATRAD  is  presented  in 
Figure  4.5.  The  various  boxes  are  self-explanatory.  More  de¬ 
tailed  flow  diagrams  of  the  transmission  function  fitting  and 
scattering  boxes  are  presented  in  Figures  4.6  and  4.7. 

In  Figure  4.7,  subroutine  SETPAR  either  sets  th* 
parameters  for  one  of  the  analytic  aerosol  size  distributions 
in  SIZDIS  or  reads  size  distribution  data  from  cards.  Sub¬ 
routine  IOR  supplies  index  of  refraction  data.  The  dashed 
box  refers  to  an  option  which  is  discussed  in  Section  4.10  and 
which  has  not  yet  been  programmed. 
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Figure  4.5  -  ATRAP  Code  Organization 
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Figure  4.6  -  Transmission  Function  Fitting 
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Figure  4.7  -  Scattering  Treatment 
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4.10  ATRAD  COMPUTATIONAL  COST;  POSSIBLE  IMPROVEMENTS 
IN  EFFICIENCY 

Two  major  realistic  test  problems,  reported  in  Chapter 
5,  have  now  been  calculated  with  ATRAD;  numerous  more  idealized 
calculations  have  also  been  completed.  From  these  calculations 
we  have  obtained  a  quite  accurate  estimate  of  how  much  computer 
time  will  be  required  for  any  particular  problem  and  of  the 
pacing  items  in  the  computation. 

The  average  amount  of  TNIVAC  1108  computer  time  re¬ 
quired  for  a  complete  run  of  ATRAD  without  aerosols  is  about 
.8  hour,  or  48  minutes.  With  cloud  aerosol  included,  the  com¬ 
puting  time  may  approximately  double,  although  for  haze 
aerosol  the  increase  may  be  smaller.  Therefore  we  are  devoting 
considerable  thought  to  improving  the  efficiency  of  the  aerosol 
computation. 

In  spite  of  any  improvements  which  we  may  make  in  the 

calculation  of  the  Mie  series  (4.26)- (4.29) ,  there  is  a 

certain  irreducible  amount  of  computation  involved  in  forming 

the  products  a  tt  ,  b  t  ,  at,  anc^  tt  •  These  four  multi- 
nn  nn  nn  nn 

plications  must  be  performed  for  every  n,  for  every  size 
parameter  a,  for  every  angle  at  which  F^  M  is  tabulated, 
for  every  level  at  which  aerosol  is  present,  and  for  every 
frequency  group.  Making  generous  estimates  of  each  of  these 
factors,  we  obtain  5  x  10 8  multiplications  for  a  complete 
problem.  At  roughly  5  ysec  per  multiplication,  this  amounts 
to  2,500  sec  or  about  .7  hours. 

We  have  two  ideas  on  reducing  the  computing  time  in 
the  Mie  calculation.  The  first  is  simply  to  use  a  much 
coarser  frequency  group  structure  for  the  Mie  calculation  than 
for  the  problem  itself,  with  some  sort  of  extrapolation  or 
interpolation  to  obtain  Mie  quantities  on  the  finer  mesh. 
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The  frequency  group  structure  of  the  problem  is  determined  by 
the  necessity,  to  resolve  rapidly  varying  quantities  such  as 
the  solar  flux  S^,  Rayleigh  scattering  coefficient 
Planck  function  B  ,  etc.  The  Mie  quantities  in  general  vary 
much  more  slowly  with  frequency. 

The  second  approach  is  to  calculate  large  tables  of 
Mie  quantities ,  store  them  on  disc  or  fast  drum,  and  then 
merely  interpolate  in  these  tables  in  subsequent  ATRAD  runs. 

New  tabular  entries  could  be  made  as  new  index  of  refraction 
data,  size  distribution,  data,  etc.,  became  available.  This 
approach  will  probably  be  implemented  in  the  near  future. 

A  further  pacing  item  in  ATRAD  is  the  matrix  inversion 
subroutine  used  for  doubling  and  for  the  foward  pass  of  the 

[81] 

Grant-Hunt  algorithm.  This  routine  was  obtained  from  Bell  Labs 
and  is  unquestionably  the.  best  general  matrix  inverter  available 
today.  However,  the  matrices  we  are  inverting  have  a  special 
form  which  is  not  being  exploited.  We  intend  to  consult  with 
a  few  of  the  leading  experts  in  linear  algebra  to  determine  if 
improvement  can  be  achieved  in  the  inversion  of  this  particular 
matrix. 
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5.  RADIATION  IN  THE  EARTH'S  ATMOSPHERE: 
ATRAD  APPLICATIONS 


This  chapter  presents  the  results  of  comparisons  be¬ 
tween  the  ATRAD  code  and  the  radiation  subroutine  of  the 
Mintz-Arakawa  global  circulation  code.  The  Mintz-Arakawa 
treatment  is  discussed  in  Section  5.1.  The  test  problem  re¬ 
sults  are  presented  in  Section  5.2. 

5.1  THE  MINTZ-ARAKAWA  RADIATION  SUBROUTINE 

The  calculation  of  radiation  heating  and  cooling  in 
the  Mintz-Arakawa  (M/A)  global  circulation  model  (GCM)  is 
necessarily  a  considerably  simplified  one  in  order  to  keep 
the  computational  work  within  practicality.  In  order  to 
determine  the  size  of  the  errors  which  result  from  this  sim- 

radiative  treatment,  a  series  of  comparisons  of 
fluxes  and  heating  rates  calculated  with  the  first-principles 
radiative  transfer  code,  ATRAD,  and  those  from  the  Mintz- 
Arakawa  two- level  GCM  are  being  made.  The  results  of  the 
first  two  comparison  cases  corresponding  to  cloud-free  con¬ 
ditions  over  land  are  reported  in  Section  5.2.  In  order  to 
make  these  comparisons  we  performed  calculations  of  radiative 
quantities  with  the  M/A  radiative  subroutine  corresponding 
to  a  particular  epoch  in  a  standard  global  circulation  calcu¬ 
lation.  (We  are  indebted  to  A.  B.  Kahle  of  RAND  for  the 
M/A  computer  code  and  for  input  data  for  it.)  In  the 
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following  section  we  outline  the  features  of  the  M/A  radiative 
model  and  describe  how  the  input  data  for  it  are  transcribed 
to  the  ATRAD  code . 

The  M/A  treatment  as  employed  in  the  climate  dynamics 
program  of  the  RAND  Corporation  is  described  by  Gates, 
et.al. ^  ,  who  give  a  definitive  description  of  the  formula¬ 

tion  and  code.  We  outline  the  principal  features  of  the  code 
here;  however,  we  are  particularly  interested  in  specifying 
how  the  M/A  atmosphere  is  constituted  in  order  that  we  may 
provide  a  proper  initialization  for  ATEAD. 

The  M/A  radiative  subroutine  determines  the  radiative 
heating  rate  in  two  layers  of  the  atmosphere  and  the  radiative 
flux  into  the  surface.  Each  of  the  layers  contains  one-half 
of  the  mass  of  the  atmosphere  lying  below  the  200  mb  level. 

The  part  of  the  atmosphere  above  200  mb  is  not  explicitly 
contained  in  the  M/A  calculation  but  is  included  in  ATRAD. 

The  relevant  radiative  quantities  to  determine  atmospheric 
heating  are  completely  specified  by  giving  the  net  radiative 
fluxes  at  the  boundaries  of  these  levels;  in  terms  of  the 
variable  a  denoting  the  fraction  of  the  atmosphere  below 
200  mb  these  are  F(a  =  0)  at  the  200  mb  level,  F(a  =  0.5) 
in  the  middle  of  the  atmosphere,  and  F (a  =  1.0)  at  the  lower 
boundary  of  the  atmosphere.  These  three  altitudes  are  in¬ 
cluded  among  the  zone  boundaries  of  ATRAD  for  which  radia¬ 
tive  fluxes  are  calculated. 

The  radiative  model  of  the  M/A  code  distinguishes 
three  spectral  regions  in  which  different  radiative  calcula¬ 
tions  are  performed.  Taking  advantage  of  the  fact  that  there 
is  little  overlap  between  the  solar  and  thermal  radiative 
sources , . the  major  frequency  division  is  between  long— wave 
and  short-wave  radiation.  The  frequency  dividing  these  re¬ 
gions  is  not  specified  in  the  M/A  code  but  is  chosen  to  be  at 
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4p  for  the  purposes  of  editing  data  from  ATRAD.  In  the  long¬ 
wave  region  the  M/A  code  takes  account  only  of  thermal  emis¬ 
sion  and  absorption  by  the  atmosphere  and  the  ground.  In  the 
short-wave  case  where  thermal  emission  is  neglected  the  spec¬ 
trum  is  divided  at  0.9u  into  two  parts;  in  .he  near  IR  region 
the  solar  radiation  is  assumed  to  be  absorbed  but  unscattered, 
while  in  the  remaining  spectral  interval  it  is  assumed  to  be 
unabsorbed  but  Rayleigh  scattered.  Using  these  spectral  in¬ 
tervals,  the  results  of  the  M/A  calculation  can  be  presented 
as  a  matrix  of  nine  net  fluxes,  corresponding  to  the  three 
spectral  intervals  and  the  three  altitudes. 

The  most  important  atmospheric  absorber  for  typical 
atmospheric  conditions  is  water  vapor.  The  M/A  calculation 
takes  account  principally  of  the  effect  of  water  vapor,  in¬ 
corporating  empirical  corrections  for  the  effect  of  C02. 

Based  on  a  single  value  q3  of  the  water  vapor  mixing  ratio 
0  ~  3/4,  and  an  assumed  saturated  vapor  pressure  above 
120  mb  in  the  stratosphere,  the  code  interpolates  a  profile 
of  water  vapor  density.  The  vertical  profile,  Q,  as  a  func¬ 
tion  of  pressure,  p,  in  mb  is  given  by 


p  ^  120  mb  , 

(5.1) 

p  <  120  mb 

P3  is  the  pressure  at  a  =  3/4  and  the  constant  k  is  deter¬ 
mined  by  matching  the  values  of  Q  at  120  mb: 
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k  = 


ln(Q3/1.7188  x  lO"6) 
ln(p3/120  mb) 


Using  this  profile  the  effective  water  vapor  amount, 
including  a  factor  for  the  pressure  broadening  of  the  absorp¬ 
tion  lines,  is  calculated  between  altitude  levels.  This 
quantity  is  subsequently  used  to  form  the  absorptivity  of  each 
layer  for  the  near  IR  short-wave  radiation  and  the  transmis¬ 
sivity  of  each  layer  for  the  IR  radiation. 

Equation  (5.1)  is  also  used  to  form  the  water  vapor 
density  input  data  for  the  ATRAD  code.  In  addition,  it  is 
necessary  to  interpolate  and  extrapolate  the  temperature. 

The  M/A  code  makes  use  of  the  calculated  values  of  temperature 
T1  in  the  upper  layer  and  T3  in  the  lower  layer.  From  these 
values  the  temperatures  at  the  200  mb  level  and  at  the  boundary 
between  the  upper  and  lower  layer  are  found  by  invoking  the 
following  functional  dependence  of  temperature  on  pressure : 


T 


<PK) 


Q  < 

6lp3 


P' 


03P1 

Pi 


where 


k  =  R/C  and 
P 


0  =  T/pK 


(5.2) 


is  proportional  to  potential  temperature.  Equation  (5.2)  in¬ 
sures  that  the  prescribed  values  of  temperature  at  p^ 

(a  =  1/4)  and  p3  (a  =  3/4)  are  recovered  and  in  addition, 
the  temperature  becomes  zero  at  p  =  0  .  In  addition,  the  air 
tempei  ature  near  the  ground  at  the  top  of  the  boundary  layer, 
T4 ,  is  also  calculated  in  the  M/A  code  from  a  consideration  of 
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the  fluxes  of  static  energy  between  the  ground  and  level  3 
(o  =  3/4) .  Subsequent  adjustments  of  the  original  input  tem¬ 
peratures  may  take  place  from  convection.  In  order  to  specify 
the  initial  temperature  profile  for  the  ATPAD  code,  we  use 
Eq.  (5.2)  from  a  =  0  to  a  =  3/4  .  To  determine  temperatures 
from  a  =  0.75  to  a  =  1.0,  we  use  the  same  expression  ad¬ 
justed  to  match  T~  and  T.: 

J  4 


T  = 


0.  -  0. 

4 _ . 

K  K 

Pa  ~  P- 


(PK) 


A  r,K 

63p4 


A  K 

04p3 


Pa  ~  P- 


P3  £  P  £  P4  •  (5.3) 


Finally,  in  the  region  above  200  mb,  where  the  M/A  code  gives 
no  guidance,  we  provide  a  rough  simulation  of  the  stratospheric 
temperature  rise,  using  an  expression  which  is  continuous 
with  Tq  at  Pq  =  200  mb  . 

(P0  +  M 

T  "  To  7F  +  kx)  '  0  £  p  £  p0  •  (5.4) 


Having  specified  the  relation  between  temperature  and 
pressure  throughout  the  atmosphere,  we  can  proceed  to  the  deter¬ 
mination  of  the  altitude  by  using  the  hydrostatic  relationship, 


where  Z4  is  the  altitude  of  the  surface  and  we  have  neglected 
the  effect  of  water  vapor  in  the  atmosphere.  In  view  of  the 
several  analytic  forms  given  above  for  the  temperature-pressure 
relationship,  we  have  approximated  it  by  linear  expressions, 

T  =  A  +  Bp  ,  within  each  of  the  ATRAD  zones  which  are  of  the 
order  of  25  mb  thick.  The  altitude  of  zone  I  becomes 
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Z(I)  -  Z(I-l)  +  | [a  In  ^TT“  +'  ~  P(I))] 


(5.5) 


where 


A  =  T(I-D  P(I)  ~  T (I )  P(I-l) 
p(D  -  P(I-l) 


and 


B  -  Td)  ~  T(I-l) 

P(l)  -  P(I-l)  * 

We  have  chosen  Z(l)  =  0,  thereby  measuring  the  altitude  rela¬ 
tive  to  the  surface. 

5 . 2  TEST  PROBLEMS 

Two  different  locations  on  the  Earth's  surface  were 
chosen  for  comparison  calculations.  The  criteria  for  choosing 
these  locations  were  first,  that  they  be  the  centers  of  grid 
squares  in  the  M/A  code,  second,  that  the  M/A  code  predict 
no  clouds  for  the  chosen  points,  and  third,  that  the  chosen 
points  have  substantially  different  water  vapor  profiles.. 

Two  such  locations  were  found  in  the  M/A  run  that  was  furnished 
to  us.  The  first  (Problem  1)  was  at  14°N,  20°E,  in  Chad  (in 
the  Sahara  desert) .  The  second  (Problem  2)  was  at  18°S,  65°W, 
in  Bolivia.  Temperature  profiles  and  water  vapor  profiles 
for  the  two  problems,  as  extracted  from  M/A  by  the  methods  of 
the  previous  section,  are  shown  in  Figures  5.1  and  5.2,  re¬ 
spectively.  Since  M/A  does  not  calculate  ozone  profiles,  an 
analytic  ozone  profile  taken  from  Green,  182  ^ 
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5.1  -  Test  Problem  Temperature  Profiles 


SSS-R-72-1255 


I 


(Z) 


We 


(z-zo)/h 


h[] 


1  +  e 


(z-zj/hy 


was  used  with  parameters  W  =  0.218  atm-cm  (total  ozone), 
h  =  4.5  km  (ozone  layer  half-width),  and  ZQ  =  23.25  km  (ozone 
maximum).  This  profile  is  shown  in  Figure  5.3. 

The  cosine  of  the  solar  zenith  angle,  ground  tempera¬ 
ture,  T  ,  and  surface  albedo  for  each  problem  were  also  taken 
directly  from  the  M/h  code.  For  Problem  1  they  are 


and  for  Problem  2, 


0.84641818 
343 . 26°K 

0  v  <  2500  cm-1 

0.2  v  >  2500  cm-1 

0.454845 
304 . 19°K 

0  v  <  2500  cm"1 

0.09  v  >  2500  cm” 1 


Note  from  Figure  5.1  that  in  both  problems  there  is  a  tempera¬ 
ture  slip  at  the  ground,  of  magnitude  31.5°K  for  Problem  1, 
and  of  magnitude  15°K  for  Problem  2. 

Both  problems  were  run  with  39  zones,  extending  from 
the  surface  to  2  mb.  A  total  of  116  frequency  groups  were 


119 


SSS-R-72-1255 


used,  32  in  the  IR  (60-2500  cm-1)  29  in  the  near  IR  (2500- 
11,000  cm-1)  and  55  in  the  "visible"  (11,000-50,000  cm"1).  A 
variable  number  of  Gaussian  angles  was  used:  6  between  60  and 
800  cm-1,  10  between  800  and  1400  cm-1,  8  between  1400  and 
2500  cm  ,  and  10  between  2500  and  50,000  cm” * .  Fewer  angles 
were  used  in  the  strong  water  vapor  and  CC>2  absorption  bands, 
based  on  sample  calculations  showing  a  high  degree  of  hemi¬ 
spherical  isotropy  in  the  intensity  in  these  bands.  Test  cal¬ 
culations  were  run  with  more  angles  in  selected  frequency 
groups  to  determine  if  the  number  of  angles  being  used  was 
sufficient;  in  no  case  did  the  fluxes  vary  by  more  than  a 
percent  upon  using  more  angles. 

Test  calculations  were  also  performed  to  ascertain 
if  the  fluxes  were  sensitive  to  variations  in  the  exponential 
fitting  process.  As  we  pointed  out  in  the  previous  report, 
the  exponents  and  coefficients  resulting  from  this  process 
are  highly  sensitive  to  perturbations  in  the  transmission 
function  data,  but  the  conjecture  was  made  that  the  fluxes 
and  intensities  would  be  insensitive  to  variations  in  the  fit. 

In  the  test  calculations,  various  small  perturbations  in  the 
transmission  function  data,  which  produced  large  perturbations 
in  the  exponential  fits,  nevertheless  produced  only  small, 
perturbations  (always  less  than  1  percent)  in  the  fluxes. 

Hence ,  our  previous  conjecture  has  been  borne  out  remarkably 
well. 

The  ATRAD  frequency-integrated  radiative  flux  F(z) 
for  Problem  1  is  given  in  Figure  5.4  as  a  function  of  pressure 
p(z) .  (NOTE:  positive  fluxes  are  downward  directed,  negative 
fluxes  are  upward  directed.)  The  Problem  1  atmosphere  is  very 
dry,  which  is  manifested  in  the  near  constancy  of  F  throughout 
most  of  the  atmosphere.  The  large  change  in  F  near  the  sur¬ 
face  can  be  explained  in  terms  of  the  large  temperature  slip 
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Figure  5.4  -  Total  (frequency-integrated)  radiative 
flux  profile  from  ATRAD,  Problem  1 
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(31.5  )  at  the  surface;  this  slip  produces  a  large  upward 
blackbody  flux  which  decreases  the  downward  flux  and  which 
penetrates  several  hundred  meters  upward  before  absorption 

and  re-emission  by  water  vapor  can  produce  a  cancelling  down¬ 
ward  flux. 

The  Problem  1  flax  spectrum  F^(z)  for  z  -  0  (top 
of  the  atmosphere)  and  z  =  zq  (surface)  is  plotted  as  a 
function  of  wavenumber  v  in  Figure  5.5.  Actually,  because 
of  the  logarithmic  wavenumber  scale,  vF^  rather  than  F  is 
plotted;  the  relationship 


/ 


vF  d (In  v) 


-/> 


F  dv 


then  insures  that  the  area  between  the  curve  and  the  horizontal 
axis  in  any  wavenumber  interval  is  really  the  integrated  flux. 
Prominent  absorption  features  due  to  CC>2  around  700  cm-1  (15y 
band),  C>3  around  1000  cm"1  (9.6y  band),  H20  in  the  nea-  IR 
(5  strong  bands  between  2500  and  12,000  cm"1),  and  03  in  the 
ultraviolet  (35,000-50,000  cm"1)  are  all  visible  in  this 
spectrum. 

Comparisons  between  ATRAD's  and  Mintz-Arakawa ' s  radia¬ 
tive  fluxes  and  "heating  rates"  are  presented  in  Table  5.1  for 
Problem  1.  The  "heating  rate"  table  simply  contains  flux 
ferences  from  the  flux  table  above  it.  The  error  in  the 
M/A  net  heating  rate  is  quite  large  for  this  problem  -  it  is 
too  small  by  about  a  factor  of  2  for  the  lower  level,  and  too 
large  by  a  factor  of  4  and  of  the  wrong  sign  for  the  upper 
level.  lhe  M/A  flux  into  the  ground  is  too  large  by  almost 
100  watts/m2,  or  about  31  percent.  The  M/A  fluxes  are  in 
general  considerably  more  accurate  than  the  M/A  heating  rates. 
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This  illustrates  the  fact  that  crude  models  to  obtain  radia¬ 
tive  fluxes  perform  reasonably  well  in  practice  but  flux 
differences  obtained  from  such  models  can  be  grossly  in  error. 
It  is  necessary  to  have  quite  an  accurate  flux  computation  if 
flux  differences  are  desired,  for  in  the  process  of  differenc¬ 
ing  one  or  even  two  significant  figures  of  accuracy  may  be 
lost. 

The  ATRAD  frequency-integrated  flux  F(z)  for  Problem  2 
is  given  in  Figure  5.6.  The  Problem  2  atmosphere  has  con¬ 
siderably  more  water  vapor  at  all  levels  than  the  Problem  1 
atmosphere  —  a  factor  of  12  more  in  total  water  content  and  a 
factor  of  19  more  at  the  surface  —  which  explains  the  slow  de¬ 
crease  of  F  with  altitude  (the  higher  in  the  atmosphere  one 
is,  the  stronger  the  upward  flux  due  to  water  vapor  emission 
becomes,  which  cancels  more  of  the  downward  solar  flux).  As 
in  Problem  1,  the  sharp  flux  change  near  the  surface  is  as¬ 
cribed  to  the  temperature  slip  (15°) .  The  change  is  much 
more  abrupt  in  Figure  5.6  than  in  Figure  5.4  because  the  air 
near  the  surface  is  so  much  wetter,  causing  the  water  vapor 
emission  to  cancel  the  temperature  slip  flux  a  very  short 
distance  above  the  surface. 

The  Problem  2  flux  spectrum  F^(z)  is  given  in  Figure 
5.7  for  z  =  0  and  z  =  zq.  For  the  same  reasons  as  in  Figure 
5.5,  vF^  rather  than  F^  is  plotted.  Note  that  the  near  IR 
H20  absorption  features  are  enhanced,  and  that  the  C02  15y 
and  9.6y  features  are  more  washed  out,  with  respect  to 
Problem  1.  This  is  because  of  the  increased  amount  of  water 
vapor  in  Problem  2. 

The.  comparison  of  ATRAD  with  M/A  for  Problem  2  is  pre¬ 
sented  in  Table  5.2.  The  M/A  and  ATRAD  fluxes  agree  much  more 
closely  for  Problem  2  than  they  did  for  Problem  1.  The  M/A 
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net  flux  errors  are  1.6  percent  for  a  =  0,  2.3  percent  for 
a  =  0.5,  and  4  percent  for  o  =  1.0.  Again,  however,  because 
of  loss  of  significance  in  subtraction,  the  heating  rate 
errors  are  substantially  larger.  The  M/A  heating  rate  in  the 
lower  level  is  too  large  by  more  than  a  factor  of  2,  and  in 
the  upper  level  is  in  error  by  7.7  percent.  The  M/A  surface 
flux  (a  =  1) ,  which  determines  the  heating  of  the  around,  is 
in  error  by  4  percent.  Thus,  the  M/A  radiation  model  makes 
very  good  predictions  of  fluxes  and  heating  of  the  lower 
level.  It  is  not  surprising  that  the  M/A  radiation  model 
performs  substantially  better  on  Problem  2  than  it  did  on 
Problem  1,  for  it  devotes  a  great  deal  of  attention  to  water 
vapor  absorption,  and  Problem  2  is  much  more  dominated  by 
water  vapor  absorption  than  Problem  1. 

It  is  clear  that  the  M/A  radiation  model  incurs  sub¬ 
stantial  error  in  the  case  of  dry  atmospheres.  However,  no 
definitive  conclusions  with  respect  to  the  usefulness  of  the 
M/A  radiation  model  for  global  climate  calculations  can  be 
drawn  from  just  these  two  problems.  Many  more  test  calcula¬ 
tions,  including  cases  with  clouds,  cases  with  realistic 
aerosol  distributions,  cases  with  high  albedo  (such  as  in  the 
Arctic),  and  nighttime  cases,  need  to  be  performed  in  order 
to  obtain  a  more  complete  picture. 

Even  when  this  picture  is  obtained,  other  factors 
need  to  be  considered.  For  example,  the  magnitude  of  the 
error  in  the  total  heating  rate,  including  latent  heating, 
adiabatic  heating,  etc.,  produced  by  errors  in  the  radiative 
contribution  must  be  ascertained.  Also,  it  might  be  the  case 
that  a  radiation  model  whose  averaged  (zonally,  diurnally, 
etc.)  predictions  are  correct  is  sufficient  for  climate  dynam¬ 
ics  calculations,  in  that  case,  it  would  be  the  average  of 
many  ATRAD  calculations  which  would  have  to  be  compared  with 
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the  average  of  many  M/A  radiation  calculations.  The  questions 
raised  in  this  paragraph  will  have  to  be  addressed  before  the 
amount  of  improvement  necessary  in  the  M/A  radiation  model 
can  be  ascertained. 
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6,  FUTURE  CONSIDERATIONS 

6.1  OROGRAPHIC  EFFECTS  ON  VERTICAL  MOMENTUM  FLUX  IN  THE 

ATMOSPHERE 

The  main  emphasis  of  the  next  contract  period  will  be 
the  development  of  three  dimensional  transient  and  steady 
state  codes.  The  parameterization  effort  will  also  continue 
attempting  to  understand  the  relationship  between  the  steady 
state  and  transient  results  and  in  arriving  at  simplified 
expressions  capable  of  being  used  in  the  global  circulation 
model. 

One  important  question  which  will  be  addressed  during 
the  parameterization  studies  is  the  desirability  of  using  an 
expression  derived  from  either  the  steady  state  or  transient 
results  as  neither  properly  expresses  the  results  truly 
representative  of  the  GCM  calculations.  The  changing  wind 
patterns  rule  out  to  a  cercain  degree  the  steady  state  re¬ 
sults  and  the  transient  results  calculated  to  this  date  have 
not  been  carried  out  late  enough  in  time  to  determine  if  they 
are  respentative  of  the  situation.  Certainly,  it  seems  worth 
while  to  pursue  this  question  on  at  least  a  semi-quantitative 
basis . 


6.2  RADIATIVE  HEATING  IN  THE  ATMOSPHERE 

In  order  to  perform  a  large  number  of  comparison  cal¬ 
culations  between  ATRAD  and  the  M/A  radiation  model,  especially 
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with  the  Mie  part  of  ATRAD  turned  on  (as  it  was  not  for  the 
problems  of  Chapter  5) ,  substantial  effort  will  have  to  be 
devoted  to  improving  the  running  time  of  ATRAD.  Some  of  our 
ideas  for  doing  so  were  discussed  in  Section  4.10.  The  main 
methods  by  which  we  intend  to  improve  ATRAD' s  running  time 
are : 

(1)  replace  the  Mie  calculation  with  tabular 
interpolation  in  large  tables  of  Mie  quan¬ 
tities; 

(2)  replace  the  exponential  fitting  calculation 
with  tables  of  coefficients  a^^  and  exponents 

V 

(3)  look  for  a  faster  matrix  inversion  subroutine; 

(4)  code  some  frequently  performed  calculations 
in  assembly  language  for  greater  efficiency. 

Once  these  improvements  are  completed,  we  intend  to 
perform  many  more  comparisons  between  ATRAD  and  M/A,  sufficient 
to  obtain  a  spatially  and  temporally  complete  picture  of  the 
accuracy  of  the  M/A  radiation  model.  Concurrently  with  this 
effort,  methods  will  be  sought  to  improve  the  M/A  code, 
especially  the  treatment  of  absorbers  other  than  H20  and  its 
crude  frequency  resolution.  Comparisons  will  also  be  made 
with  the  radiation  subroutines  of  other  Global  Circulation 
Codes  (it  is  likely  that  the  GCM  of  the  National  Center  for 
Atmospheric  Research  will  be  first) .  These  calculations  will 
provide  a  basis  for  intercomparison  of  different  GCM  radia¬ 
tion  subroutines.  Subsequently ,  a  more  comprehensive  assess¬ 
ment  of  radiative  subroutine  accuracy  and  efficiency  can  be 
made. 
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APPENCiX  A  -  ATRAD  INPUT 

The  input  for  ATRAD  consists  entirely  of  FORTRAN 
NAMELIST's,  which  are  the  most  convenient  and  versatile  way 
of  performing  data  input  in  the  FORTRAN  language.  NAMELIST's 
also  have  the  advantage  of  being  exceptionally  easy  for  even 
non-programmers  to  use,  so  that  the  code  could  be  exercised 
by  personnel  with  a  minimal  amount  of  instruction. 

There  are  two  possible  deck  set-ups  for  running  ATRAD. 
The  first  (the  "normal"  mode)  is  used  for  setting-up  and 
starting  a  problem  ab  initio,  and  is  illustrated  in  Figure 
A.l.  The  second  (the  "restart"  mode)  is  used  for  restarting 
a  problem  from  a  dump  tape  at  a  predetermined  frequency  group, 
and  is  illustrated  in  Figure  A. 2. 

The  variables  required  by  each  NAMELIST  are  described 
below.  If  the  variable  is  not  specified  in  the  NAMELIST,  it 
will  assume  a  default  value  which  is  given  in  parentheses  at 
the  end  of  the  description  of  the  quantity.  The  subscript 
convention  employed  is  as  follows: 

I  is  a  level  index, 

J  is  an  angular  index,  and 
NU  is  a  frequency  index. 
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Figure 


WISCOM  NAM2LIST 

Input  Control  Parameters 
(listed  in  Section  A.l). 

END 


Only  if 
PRSETZ  =  TRUE 


FREQS  NAMELIST 

WAVNUM,  NNU 
(listed  in  A.  3) 

END 


\  Only  if 

PRSETF  =  TRUE 


LEVELS  NAMELIST 

Z,  NZ,  IBD,  NBD 
(listed  in  A.  2) 

END 


STRUCT  NAMELIST 

Pf  T,  H 2 ODEN,  03DEN, 
AERDEN,  NAER,  NMAT 
(listed  in  A.  4) 


Only  if  NOPT  =  0 


END 


AEROS  NAMELIST 

RAD,  AERNUM,  NDAT 
(listed  in  A. 5) 

END 


J  Only  if  NAER (I)  =  0 
(  or  -1  for  at  least 
’  one  I . 


A.l  Input  Deck  Set-up  for  a  Normal  Run  of  ATRAD. 
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RSTART  NAMELIST 


NUSTRT,  INPT 
(listed  in  A. 6) 


END 


WISCOM  NAMELIST 


Input  Data 
(listed  in  A.l) 


Only  if  INPT  >  0 


END 


AEROS  NAMELIST 


RAD,  AERNUM,  NDAT 
(listed  in  A. 5) 


END 


Only  if  NUSTRT  «  1 
and  NAER(I)  =  0  or 
-1  for  at  least- 
one  I. 


Figure  A. 2  Input  Deck  Set-up  for  Restarting 
ATRAD  from  Dump  Tape  at  Frequency 
Group  NUSTRT. 
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A.l  WISCOM  NAMELIST 


Variable  Type 
PR1  Logical 

PR2  Logical 


PR3  Logical 


PR4  Logical 


SCTPRT  Logical 


Description 

If  TRUE,  atmospheric  structure  (P, 

T,  etc.)  is  edited.  (FALSE) 

If  TRUE,  edits  frequency-integrated 
fluxes,  heating  rates,  and  intensi¬ 
ties,  as  well  as  the  flux  spectrum 
at  the  upper  and  lower  boundaries 
of  the  atmosphere.  (TRUE) 

If  TRUE,  edits  fluxes,  heating  rates, 
intensities,  incident  solar  flux, 
and  earth-atmosphere  albedo  for  each 
frequency  group.  (FALSE) 

If  TRUE,  edits  solar  flux,  Planck 
function,  scattering  coefficient, 
and  continuum  absorption  coefficient 
for  both  the  current  and  previous 
frequency  group,  and  the  percent 
change  of  each  quantity  from  the 
previous  frequency  group.  Also  edits 
the  continuum  part  of  the  optical 
depth  of  each  zone,  and  various  sur¬ 
face  quantities  such  as  the  direction¬ 
al  emissivity  e(y)  and  reflection 
matrix  rG(y,y').  (FALSE) 

If  TRUE,  edits  truncation  informa¬ 
tion,  Mie  scattering  and  absorption 
coefficients,  Mie  phase  function 


r* 
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Variable 

SCTPRT 
(cont. ) 

MIEPR(l) 


MIEPR (2) 


MIEPR (3) 


gype 

Logical 


Logical 


Logical 


Logical 


Description 

both  before  and  after  renormaliza¬ 
tion,  Rayleigh  and  total  scattering 
coefficients,  and  total  phase  func¬ 
tion,  (FALSE) 

If  TRUE,  euivs  on  each  call  to 
MIESCT  the  real  and  imaginary  parts 
°f  An,  an,  and  bn,  which  are  coeffi¬ 
cients  in  the  Mia  series.  This 
edit  is  primarily  useful  when  the  Mie 
package  is  being  driven  separately; 
it  should  never  be  turned  on  for  a 
normal  radiation  calculation,  as  it 
would  produce  a  vast  amount  of  print¬ 
ing.  (FALSE) 

If  TRUE,  edits  from  each  call  to 
MIESCT  the  Mie  phase  function  and  the 
efficiency  factors  and  cross-sections 
for  scattering,  absorption,  and 
extinction.  Same  comments  apply  as 
for  MIEPR (1).  (FALSE) 

If  TRUE,  edits  information  relevant 
to  integration  over  aerosol  size 
distribution  in  MIE  and  prints  Mie 
phase  function  and  corss-sections 
for  absorption  and  scattering  re¬ 
sulting  from  that  integration. 

(FALSE) 
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Variable  Type 
FITPRT  Logical 


FITPR2  Logical 


TPR1  Logical 


TPR2  Logical 


Description 

If  TRUE,  edits  for  each  frequency 
group  and  for  each  molecular  species 
(H2°'  ^°2+'  ^3)  t^le  transmission 
function,  its  exponential— sum  approx— 
imant,  and  the  percent  error.  Also 
edits  the  coefficients  a^  and  expo¬ 
nents  of  the  approximant. 

(FALSE) 

If  TRUE,  edits  the  exponential-sum 
fitting,  including  the  coefficients 
ai  and  exponential  factors  0^  the 
residual,  and  the  polynomial  (4.38), 
for  each  iteration.  Also  edits 
information  about  pairs  of  exponen¬ 
tial  factors  which  are  combined. 
(FALSE) 

If  TRUE,  performs  a  "short"  edit  for 
each  call  to  TRNPRT  (Grant-Hunt 
algorithm) ,  including  the  single- 
scattering  albedo,  line  absorption 
coefficient,  optical  depth,  doubling 
parameters,  source  vectors,  and 
resultant  (diffuse)  intensities  for 
each  zone.  (FALSE) 

If  TRUE,  performs  all  of  the  edits 
of  TPR1  plus  the  Planck  and  solar 
source  vectors,  the  reflection  and 
transmission  matrices  for  each  zone, 
and  the  various  vectors  and  matrices 
(E,  etc.)  used  in  the  Grant-Hunt 
algorithm.  (FALSE)  145 
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Variable 

PRSETZ 


OZONE 


TROPO 


NCLOUD 

CLDBAS (1,2) 

CLDTOP (1,2) 


Type 

Logical 


Logical 


Real 


Integer 


Description 

If  TRUE,  the  user  intends  to  input 
the  zone  structure  through  the 
"LEVELS"  NAMELIST.  If  FALSE,  user 
intends  to  set  OZONE,  TROPO,  NCLOUD, 
CLDBAS,  CLDTOP,  NZONES,  and  EXPAN 
(see  below)  with  which  the  code  will 
calculate  the  zone  structure. 

(FALSE) 

If  TRUE,  place  the  uppermost  level 
at  50  km  altitude.  If  FALSE,  take 
the  uppermost  level  at  TROPO  km 
(the  tropopause  or  some  level  just 
below  the  ozone  layer) .  In  the 
latter  case,  the  incident  solar  flux 
will  be  truncated  at  03CUT.  (TRUE) 

If  OZONE  =  FALSE,  the  height  in  km 
of  the  uppermost  level  in  the  prob¬ 
lem  is  set  to  TROPO.  This  level 
should  be  reasonably  close  to  the 
bottom  of  the  ozone  layer.  (15. 0) 

The  number  of  cloud  layers  (only  0, 
1,  or  2  are  allowed) .  (0) 


Real  Heights  in  kin  of  bases  of  lower  (1) 

and  upper  (2)  cloud.  (2*0.) 

Real  Heights  in  km  of  tops  of  lower  (1) 

and  upper  (2)  cloud.  (2*0.) 


r'  146 

A- 7 


SSS-R-72-1255 


Variable 
NZONES (1-5) 


EXPAN (1-5) 


PRSETF 


NWAV 

WAV (1-15) 


Type 

Real 


Real 


Logical 


Integer 

Integer 


Description 

Number  of  zones  in  each  of  the  sub- 
regions  into  which  the  clouds  parti¬ 
tion  the  atmosphere.  E.g.,  NZONES(l) 
is  the  number  of  zones  from  the 
surface  to  the  first  cloud  base,  or, 
if  NCLOUD  =  0-,  is  the  total  number 
of  zones.  NZONES  (2)  is  the  number 
of  zones  in  cloud  1,  etc.  (20,4*0.) 

Expansion  factors  for  each  sub- 
region.  E.g.,  in  sub-region  1,  each 
zone  has  a  geometrical  width 
EXPAN (1)  times  the  width  of  the  zone 
below  it.  This  enables  one  to  gen¬ 
erate  zones  of  more  or  less  constant 
mass,  increasing  mass,  etc. 

(1.1, 4*0.) 

If  TRUE,  user  intends  to  input  fre¬ 
quency  group  structure  through 
"FREQS"  NAMELIST.  If  FALSE,  user 
intends  to  set  NWAV,  WAV,  DWAV  or 
NGRPS,  IWV  (see  below)  and  let  the 
code  calculate  the  frequency  group 
structure.  (FALSE) 

Number  of  frequency  "regions."  (12) 

Wavenumber  boundaries  in  cm*"^  for 
each  frequency  "region."  E.g., 

WAV(l)  <  WAV (2)  are  the  wavenumbers 
bounding  frequency  region  1.  Each 
value  must  be  an  integer  multiple 
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Variable 


Description 


WAV (1-15) 
(cont . ) 


DWAV  (1-14) 


NGRPS 


IWV  (1-30) 


Type 

Integer 


Integer 


Integer 


Integer 


of  20  cm”1.  (60 ,  600,  760,  1300, 

1000,  2500,  3200,  4400,  7400,  9000, 
32000,  35000,  50000) 

Frequency  group  width  in  cm  1  for 
all  of  the  groups  within  a  frequency 
"region"  defined  above.  Each  value 
must  be  an  integer  multiple  of 
20  cm  1.  E.g.,  between  WAV(l)  and 

WAV (2)  all  frequency  groups  are  of 
sise  DWAV (1) .  (60,  40,  60,  100,  140, 

140,  200,  300,  400,  500,  1000,  1500) 

If  15  >  NGRPS  >  0,  it  overrides  the 
WAV,  DWAV  option.  The  NGRPS  option 
allows  one  to  do  NGRPS  disjoint  fre¬ 
quency  groups  with  boundaries  IWV 
(see  below) ;  it  can  be  used  to 
sample  the  spectrum.  (0) 

Used  when  NGRIS  >  0.  Successive 
pairs  of  IWV  values  define  the  fre¬ 
quency  boundaries  (cm  J")  of  the  fre¬ 
quency  groups  being  sampled.  E.g., 
if  NGRPS  =  2,  the  first  group  is 
[IWV (1) , IWV (2) ]  and  the  second  group 
is  [IWV (3) ,IWV (4) ] .  (30*0) 

Wavenumber  boundaries  in  cm  1  deli¬ 
neating  regions  in  which  the  number 
of  quadrature  angles  is  constant. 

This  allows  one  to  use  fewer  angles 
in  frequency  regions  where  the 
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Variable 

IWANG  (1-15) 
(cont. ) 

NUMANG (1-14) 


NUMANG (15) 

RADOW (1-14 ) 


NOPT 


Type 

Integer 


Integer 


Integer 


Logical 


Integer 


Description 

intensity  is  more  nearly  isotropic 
(in  each  hemisphere  separately) . 

(50,  50000,  13*0) 

Number  of  quadrature  angles  used  in 
each  hemisphere  in  the  corresponding 
frequency  region  delineated  by 
IWVANG.  E.g.,  the  number  of  quadra-' 
ture  angles  between  IWVANG (1)  and 
IWVANG (2)  is  NUMANG  (1)  .  (6,  13*0) 

Number  of  wavenumber  regions  der- 
lineated  by  IWVANG.  (1) 

If  RADOW (K)  =  TRUE,  then  use  Radau 
quadrature  in  wavenumber  interval 
[TWANG  (K) ,  IWANG  (K+l)  ]  .  Otherwise, 
Gaussian  quadrature.  (14*FALSE) 

Atmostpheric  structure  flag. 

=  0  Input  atmospheric  structure 
through  "STRUCT"  NAMELIST. 

=  6  Calculate  structure  from  user- 
supplied  analytic  forms  con¬ 
tained  in  subroutines  TEMPER, 
SPFHUM,  03D,  AERD,  NMATL, 

NAERO  (default  forms  of  these 
subroutines  are  used  if  user 
does  not  intervene.) 

=  7  Restart.  Read  in  structure 
from  dump  tape. 
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Variable 

NOPT 
(cont . ) 


NOPT1 

NOPT2 


PO 

DAY 

TIME 

LONG 

LAT 

CUTPLK 


Type 

Integer 


Integer 


Integer 


Real 

Real 

Real 

Real 

Real 

Real 


Description 

Otherwise,  use  one  of  the  following 
standard  atmospheres: 

=  1  Tropical 
=  2  Mid-latitude  summer 
=  3  Mid-latitude  winter 
=  4  Sub-arctic  summer 
=  5  Sub-arctic  winter 

Phase  function  renormalization  flag. 

=  1  Use  Grant  method. 

=  2  Use  Wiscombe  method.  (2) 

If  >0,  stop  calculation  just  before 
entering  frequency  loop.  This  al¬ 
lows  user  to  check  that  his  problem 
set-up  is  correct  before  doing  a 
full  run  of  the  code.  (0) 

Surface  pressure  in  mb  for  case 
NOPT  =  6.  (1013.) 

Day  of  the  year,  January  1  being 
day  1.  (1.) 

Greenwich  time,  in  hours.  (0). 

Longitude,  in  degrees,  counted  posi¬ 
tive  west  of  Greenv/ich.  (117.) 

Latitude,  in  degrees.  (33.) 

Wavenumber  in  cm  ^  above  which 
Planck  function  is  set  to  zero. 
(3333.) 

iso 
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Variable  Type 

CUTSOL  Real 

CUTRAY  Real 


03CUT  Real 


NMATG  Integer 

TG  Real 

WNSURF (1-11)  Real 


NSURF (1-10 )  Integer 


Description 

Wavenumber  in  cnf^  below  which 
incident  solar  flux  is  set  to  zero. 
(1666.) 

Wavelength  in  microns  above  which 
Rayleigh  scattering  is  neglected. 

(3.) 

Wavenumber  in  cm-1  above  which  inci¬ 
dent  solar  flux  is  set  to  zero  in 
case  OZONE  =  FALSE.  (33333.) 

Surface  material  flag.  (1) 

Surface  temperature  in  °K.  (300.) 

Frequency  boundaries  in  cm-1  speci¬ 
fying  regions  within  which  the 
surface  boundary  condition  is  calcu¬ 
lated  according  to  the  corresponding 
value  of  NSURF.  (50. ,50000. , 9*0. ) 

Between  WNSURF (K)  and  WNSURF (K+l) , 
the  surface  boundary  condition  is 
flagged  by  NSURF (K) .  The  options 
are: 

NSURF  =  1  Hemispherical  reflectivity 
or  hemispherical  emissivity 
supplied  by  user  as  function  of 
wavelength,  in  subroutine 
SURF1.  Diffuse  emission  and 
reflection  assumed. 

bf 
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Variable 

NSURF (1-10) 
(cont . ) 


NSURF (11) 

DELO 

NOTHG (I) 


TyPe  Description 

Integer  =  2  Directional  emissivity  or 

directional-hemispherical 
reflectivity  supplied  by  user 
as  a  function  of  angle  and 
wavelength,  in  subroutine  SURF 2 . 
Diffuse  reflection  assumed. 

=  3  Azimuthally-averaged  bidirec¬ 
tional  reflectivity  supplied  by 
user  as  a  function  of  angle  of 
incidence,  angle  of  reflection, 
and  wavelength,  in  subroutine 
SURF 3. 

(1,9*0) 

Integer  Number  of  non-zero  WNSURF  entries. 

(2) 


Real  Convergence  criterion  for  aerosol  size 

integration  in  MIE.  (0.001) 

Logical  if  NOTHG  (I)  =  TRUE,  use  the  normal 

procedure  to  calculate  the  Mie  phase 
function  for  zone  I.  If  NOTHG  (I)  = 
FALSE,  calculate  the  Mie  phase  func¬ 
tion  only  at  0°,  calculate  the  param¬ 
eter  g  such  that  the  Henyey- 
Greenstein  phase  function  matches 
that  value,  and  use  the  HG  phase 
function  at  all  other  angles. 

(40*FALSE) 
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Variable 
NSTEP ( 1-6) 


NDOUB 

DANGLO 


DANGHI 

DANGMX 


Description 

Integer  Used  in  the  determination  of  the 

angular  mesh  on  which  Mie  phase 
function  is  calculated.  A0  is 
calculated  in  MIE,  and  from  it  the 
angular  mesh  is  calculated  as  fol¬ 
lows:  NSTEP (-1)  steps  of  A0, 

NSTEP (2)  steps  of  2A0,  NSTEP (K) 

steps  of  2K-1A  0 ,  -  All  values 

of  NSTEP  must  be  even.  (16,3*10,2*0) 

Integer  Number  of  non-zero  entries  in  NSTEP. 

(4) 


Lower  boundary  on  A0  in  degrees 
(see  NSTEP) .  Must  be  a  multiple 

of  ^®tab'  the  interval  for  which 
Tn  are  tabulated.  (0.1) 

Upper  boundary  on  A 6  in  degrees. 

Must  be  a  multiple  of  A0  _  (see 

tab 

DANGLO).  (1.0) 


Real 


If  at  any  point  in  the  NSTEP  proced¬ 
ure  described  above,  2K-1A0  > 

DANGMX,  then  the  NSTEP  procedure  is 
terminated  and  the  rest  of  the  steps 
to  90°  are  taken  with  an  increment 
A@m  which  is  as  close  as  possible  to 

DANGMX,  is  a  multiple  of  A0^  .  (see 

tab 

DANGLO) ,  and  yet  allows  an  even 
number  of  steps  (as  required  by  the 
Simpson-rule  integration  used  to 
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in  MIE  and  to  truncate  it  in  TRUNCT) 
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Variable 

ANGMAX 


CUTMIN 


NPTS 


MINPTS 


THMN, 

THMX 


Type 

Real 


Real 


Integer 


Description 

The  maximum  allowed  angle,  in  degrees, 
at  which  the  Mie  phase  function  may 
be  truncated  (TRUNCT  chooses  pro¬ 
gressively  larger  angles  in  at¬ 
tempting  to  perform  a  satisfactory 
truncation).  -(20.) 

The  maximum  allowed  value  at  0°  of  the 
truncated  Mie  phase  function  (if  the 
TRUNCT  search  reaches  ANGMAX  without 
satisfying  the  CUTMIN  criterion,  the 
default  truncation  at  ANGMAX  may  pro¬ 
duce  a  P^  (o°)  >  CUTMIN).  (5.) 

The  maximum  number  of  transmission 
function  data  points  to  be  used  in 
the  exponential  fitting  procedure. 

(50) 


Integer  The  minimum  number  of  transmission 

function  data  points  to  be  used  in 
the  exponential  fitting  procedure. 
(10) 


Double  The  lower  and  upper  limits,  0  .  and 

Precision  .  , .  ,  .  raln 

“max'  on  t^ie  search  interval  used  to 

find  a  new  exponential  factor  0 
in  the  exponential-sum  fitting 
iteration.  These  limits  are  im¬ 
posed  because  0=0  and  0=1 
correspond  to  unphysical  values 
of  the  exponent  k  (°°  and  9,  re¬ 
spectively).  (l.D-8,. 99999999 DO) 
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Variable 

TRMIN 

TRMAX 


NSRCH 


NSRCH2 


COALES 


Type  Description 

Rea-1  The  smallest  value  of  the  transmission 

function  to  be  used  as  a  data  point 
for  exponential  fitting.  (.005) 

Real  If  the  transmission  through  the 

entire  vertical  mesh  on  a  slant  path 
specified  by  SLANT  is  >  TRMAX  for  a 
given  molecular  species,  then  the 
transmission  for  that  species  is  taken 
to  be  unity  and  no  fitting  is  done. 
(.99) 


Integer 


Integer 


Double 

Precision 


The  maximum  number  of  search  points 
taken  on  f6min/emax3  in  locating  a 
new  exponential  factor  6  in  the 
exponential  fitting  iteration.  (1001) 


The  number  of  search  points  taken  in 
intervals  [ -  A6,  6i  +  A6]  around 
each  of  the  old  exponential  factors 
0j,  in  the  process  of  finding  a  new 
exponential  factor  6  (since  exponen¬ 
tial  factors  often  occur  in  close 
pairs,  these  intervals  should  be  sub¬ 
jected  to  a  more  refined  search) . 

(21) 


The  criterion  for  whether  or  not  to 
coalesce  a  close  pair  of  exponential 
factors  If. 


An  0.+1  -  An  0± 


Jin  6j_+i  +  &n 


<  COALES 
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Variable 

COALES 
(cont. ) 

CTEST 


RTEST 


ITMAX 

ISEC1 

SLANT 


Double 

Precision 


Double 

Precision 


Double 

Precision 


Description 

then  e^/0i+1  are  replaced  by  a  single 
exponential  factor.  (.05) 

If  any  coefficient  a^  in  an  exponen¬ 
tial-sum  fit  is  <  CTEST,  the  corres¬ 
ponding  term  in  the  fit  is  dropped. 
This  number  should  not  be  much  larger 
than  the  estimated  error  in  the  trans¬ 
mission  function  data.  (.001) 


The  exponential  fitting  iteration  is 
stopped  when 


^hew  “  Rold 
R 


new 


<  RTEST 


where  R  is  the  residual  and  'old' 
and  'new1  refer  to  the  previous  and 
current  iterations.  (l.D-16) 


Integer  The  maximum  allowed  number  of  itera¬ 

tions  in  the  exponential-sum  fitting 
procedure.  (150) 


Integer  Maximum  number  of  iterations  of 

secant  method  used  in  pair-coalescence 
procedure  in  FITTRN .  (10) 


If  uv  is  the  vertical  absorber 
amount  of  a  given  molecular  species, 
then  the  smallest  value  of  the  trans¬ 
mission  function  T^  (u)  used  in 
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Variable  Type 

SLANT  Real 

(cont. ) 


CHGMAX  Real 


FACT  Real 


SCTMIN  Real 


Description 

fitting  is 

TAv  (SLJiNT*uv) 

provided  this  value  is  >  TRMIN.  (2.0) 

If  >  0,  overrides  the  normal  surface 
boundary  condition  flags  NSURF  and 
WNSURF  and  sets  the  albedo  =  CHGMAX 
(for  zero  albedo,  input  CHGMAX  = 
10"30).  (0.) 

The  fraction  of  the  maximum  primary 
layer  optical  depth,  ATmax,  which  is 
actually  used  as  a  doubling  interval 
in  the  Grant-Hunt  algorithm  in 
TRNPRT .  (0.5) 

If  the  single-scattering  albedo 
to  <  SCTMIN,  then  the  part  of  sub¬ 
routine  TRNPRT  which  neglects  scatter¬ 
ing  entirely  (involving  considerably 
less  calculation)  is  used.  (l.E-5) 


A. 2  LEVELS  NAMELIST 


Z (I)  Real  The  heights  above  the  surface  of  the 

various  levels  in  the  vertical  mesh, 
starting  from  the  surface  £Z (1 ) ] ,  in 
km. 

NZ  Integer  The  number  of  levels. 


rrr 

/ 
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Variable  Type 

IBD(l-5)  Integer 


NBD  Integer 


Description 

The  indices  of  the  levels  which  - 
divide  the  mesh  into  sub— regions . 
E.g.,  Z(IBD(1))  is  the  height  of  the 
lower  cloud  basef  and  Z(IBD(2))  is 
the  height  of  the  lower  cloud  top. 

If  there  are  no  clouds,  IBD(l)  =  NZ. 

The  number  of  sub-regions  required  by 
the  presence  of  clouds 
(NBD  =  2  x  #  clouds  +  1) . 


A. 3 


FREQS  NAMELIST 


WAVNUM (NU)  Integer 


NNU 


Integer 


The  frequencies  in  cm-1  bounding  the 
frequency  groups,  in  increasing  order 
(WAVNUM (1)  is  the  smallest.) 

The  number  of  non-zero  WAVNUM' s  (or 
one  plus  the  number  of  frequency 
groups . ) 


A. 4  STRUCT  NAMELIST 


P(I)  Real 

T(I)  Real 

H20DEN (I)  Real 


Pressures  in  mb  corresponding  to 
Z(I).  (P  (1 )  is  the  surface  pressure.) 

Temperatures  in  degrees  Kelvin.  T(l) 
is  the  temperature  of  the  air  immediate¬ 
ly  above  the  ground,  and  may  be  dif¬ 
ferent  from  TG,  the  ground  tempera¬ 
ture. 

Water  vapor  density  in  g/m3. 
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Variable 

Type 

Description 

03DEN  (I) 

Real 

Ozone  density  in  atm-cm/km. 

AERDEN ( I ) 

Real 

Aerosol  number  density  in 

3 

partxcles/cm  . 

NAER(I) 

Integer 

Flag  specifying  aerosol  size 

distribution: 

=  -2  Use  size  distribution  from 
level  above. 

=  “1  Card  input,  histogram  form 

=  0  Card  input,  pointwise  form 

=  1  Deirmendjian  Haze  M 

=  2  Deirmendjian  Haze  L 

=  3  Deirmendjian  Haze  H 

=  4  Deirmendjian  Cloud  C.l 

=  5  Deirmendjian  Cloud  C.2 

=  6  Deirmendjian  Cloud  C.3 

=  7-15  Various  Gaussians 

=  16-20  Various  log-normals 

NMAT(I)  Integer  Flag  specifying  aerosol  material: 

=  1  Water 

=  2  Sahara  dust  (Volz) 

=  3  Dust,  y  (Volz) 

=  4  Sea-salt  (Volz) 

=  5  Water  solubles,  B1  (Volz) 

=  6  Water  solubles,  M  (Volz) 

=  7  Water  solubles,  T2  (Volz) 

=  8  Soot 

=  9  Ice 

A-20 
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A. 5  AEROS  NAMELIST 


Variable  Type 

RAD(K,I)  Real 

NDAT(I)  Integer 

AERNUM (K, I )  Real 


Description 

Finite  mesh  of  aerosol  radii,  in 
microns,  K  =  1  to  NDAT(l),  on  which 
aerosol  size  distribution  AERNUM  is 
specified  (level  I) . 

Number  of  values  of  RAD  for  level  I. 

Aerosol  number  density  distribution 
in  particles/ cm^/micron  (level  I). 
NAER(I)  =  0:  Pointwise  data,  (RAD (K, I ) , 
AERNUM (K, I ) ) ,  K  =  1  to  NDAT(I). 
NAER(I)  =  -1:  Histogram  data,  value 
of  distribution  =  AERNUM (K, I)  between 
RAD(K,I)  and  RAD(K+1,I). 


A. 6  RSTART  NAMELIST 


NUSTRT  Integer 


INPT  Integer 


The  number  of  the  frequency  group  with 
which  the  calculation  is  to  begin  (if 
NUSTRT  =  1,  this  saves  the  expense 
of  setting  up  a  problem  again,  but 
otherwise  does  the  problem  from 
scratch) . 

If  >  0,  WISCOM  NAMELIST  follows 
immediately  containing  changes  to 
one  or  more  input  variables.  If  =  0, 
WISCOM  is  edited  only.  If  <  0,  none 
of  the  above. 
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APPENDIX  B  -  ATRAD  OUTPUT 


In  order  to  illustrate  the  output  of  ATRAD,  a  sample 
calculation  has  been  performed  for  frequency  group  29,000- 
29,800  cm  .  a  De.i  rmendjian  Hase  M  water  droplet  aerosol 
was  used  at  one  level  only  (the  surface  level) ;  no  aerosol 
was  assumed  at  any  other  levels.  Selections  from  the 
printed  output  (there  is  also  tape  output)  are  presented 
in  the  following  pages  and  described  below  in  order  of 
appearance: 

B.l  EDITS  DONE  BEFORE  LOOP  OVER  FREQUENCIES 

(a)  The  input  parameters  in  NAMELIST  WISCOM  (see 
Appendix  A,  Section  A.l) ; 

(b)  The  atmospheric  structure  (see  Appendix  A, 
Section  A. 4) ; 

(c)  The  quantities  F ,  DL,  and  UMAX,  which  are  re¬ 
spectively  the  integrands  in  the  expressions  (4.22)  as  a 
function  of  level,  the  effective  absorber  amounts  Au  for 
each  zone  (integrals  of  the  F's),  and  the  sum  of  all  the 
Au's  .  All  these  quantities  are  given  as  a  function  of  ab¬ 
sorbing  species  using  the  LOTRAN  prescriptions. 

(d)  The  solar  zenith  angle  (MUSUN)  and  earth-sun 
distance  (RSUN) ; 
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B • 2  FREQUENCY  group  29  . nnn-?Q r ann 

(e)  Tha  quadrature  angles  (MU(J))  and  weights  (WT(J)); 

(f)  Edits  from  the  integration  over  aerosol  sizes 

in  subroutine  MIE.  Successive  Romberg  results  for  SGSCT 
(proportional  to  3..  M)  and  SGEXT  (proportional  to  3  + 

av ,M  and  Pv,M(0°>  ^Pv,M  is  filled  in  at  other  angles  with 
the  Henyey-Greenstein  phase  function)  are  followed  by  an  edit 

of  Pv,M  (pHASEM)  for  all  ANGL ( 6)  and  XMU(cose),  and  by  the 

integrals  of  n(a,z)  times  each  of  o  ,  o  ,  ,  and  o 

sea  abs'  ext 

(average  cross-sections) .  Finally  the  renormalization  factor 
(4.35)  is  printed. 

(g)  Mie  phase  function  truncation  information,  in¬ 
cluding  value  of  P^M(0°)  (AX)  and  factor  F  (FTRUNK)  . 

Oi)  BETASC  (3v,m),  BETABS  (a^),  vectors  S0URC1  and 
S0URC2  (Pv(z,y^#yQ)  and  pv 2  /-yo) )  ,  and  matrices  PHFCN1 

and  PHFCN2  (pv /M ( z / /  M j )  and  pv  ,M(Z  >  •  “Kj  ) )  before  and  after 
renormalization  by  the  methods  of  Section  4.7.7. 

(i)  The  results  of  the  transmission  function  fitting 
process  for  C>3 ,  the  only  active  absorber.  Only  the  special 
case  of  a  1-term  fit  is  used  because  the  total  optical  depth 
of  0^  is  small. 

(j)  A  "short"  edit  of  the  single  pass  (because  of 
the  single  term  exponential  fit)  through  the  Grant-Hunt 
algorithm. 

(k)  The  Planck  function  B^(B),  scattering  coefficient 

3V  (BETASC),  continuum  absorption  coefficient  aCOnt  (BETABS), 

and  solar  flux  (SOLFLX)  for  the  current  and  previous  frequency 

group.  Also  the  optical  depth  DTAU(CONT)  of  each  zone  due  to 
cont 
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(l)  Boundary  condition  information:  emissivity  e' 

(EPS),  e^B^(Tg)  (W)  ,  (v*£/Vi0)uo  S^/ir  (Wl)  ,  and  matrix 

rG  =  2V> j  cj  Pv"  (Vi/Vj)  (RG)  . 

(m)  Fluxes  and  intensities  for  each  level  and  heating 
rates  for  each  zone. 

B . 3  SUMMARY  PRINTS 

The  summary  prints  contain: 

(1)  the  frequency-integrated  flux  F(z)  at 
each  level  and  the  frequency-integrated 
heating  rate  in  each  zone. 

(2)  the  flux  spectrum,  F^  ,  as  a  function  of 
frequency  group,  at  the  upper  boundary  of 
the  problem  and  at  the  surface. 
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A  DETAILED  RADIATION  MODEL  FOR  CLIMATE  STUDIES; 
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1.  INTRODUCTION 

An  atmospheric  radiative  transfer 
computer  code  lias  been  developed,  tested 
and  applied  to  the  calibration  of  the 
radiative  transfer  subroutine  of  one  of 
the  climate  dynamics  computer  codes.  The 
code  to  be  described  below  has  the  objec¬ 
tive  of  calculating  radiative  fluxes  of 
short  and  long  wave  radiation  in  ord  ;r  to 
de termine- highly  accurate  atmospheric 
heating  rates.  Through  this  formulation, 
which  takes  account  of  frequency  and  an¬ 
gular  dependence,  scattering  from  mole¬ 
cules  and  aerosols  ,  and  general  surface 
boundary  conditions,  we  hope  to  make 
available  a  radiative  transfer  benchmark 
code  serving  as  a  standard  for  comparison 
with  calculations  containing  less  compre¬ 
hensive  treatments. 

In  several  intensively  investi¬ 
gated  problems  of  atmospheric  motion  the 
role  of  radiative  transfer  is  an  important 
one.  One  of  these  is  the  global  circula¬ 
tion  of  the  atmosphere  in  which  short  and 
long  wave  radiation  form  a  major  subsys¬ 
tem  of  the  strongly  interactive  radiative- 
fluid  dynamic  atmospheric  system.  In  in¬ 
vestigations  of  climate  the  geographical 
and  seasonal  dependence  of  the  heating  of 
the  atmosphere  interacting  with  surface 
albedo  and  cloud  cover  provides  the  source 
of  energy  for  the  atmospheric  circulation. 
Another  active  y  investigated  field  of 
fluid  dynamics  in  which  radiation  plays  a 
role  is  the  planetary  boundary  layer. 
Diurnal  changes  within  this  layer  give 
rise  to  winds  and  clouds  having  profound 
influence  on  Man.  To  mention  a  few  of  the 
phenomena  induced  by  the  heating  of  the 
atmosphere  near  the  ground  there  are 
slope  winds,  land-sea  breezes,  heat  is¬ 
land  effects  and  radiation  and  advcction 
fogs.  Finally,  on  the  micro-scale  the 
environment  of  plants  and  animals  is 
largely  determined  by  the  radiation  prop¬ 
erties  of  their  immediate  surroundings. 

A  number  of  t lie  phenomena  mention¬ 
ed  above  have  been  the  subject  of  detailed 
theoretical  investigations  within  the  last 
decade.  If  the  results  of  these  studies 
are  to  be  oi  value  in  a  predictive  sense 
it  is  clear  that  the  treatment  of  radia-  j* 
five,  transfer  must  be  sufficiently  quan-  J 


titativc  to  determine  the  time  ar d  space 
dependence  of  the  atmospheric  heating  rate. 

We  shall  return  to  the  question  of 
how  accurate,  in  principle,  a  theoretical 
description  of  atmospheric  radiative 
transfer  can  be  made.  Even  within  the 
framework  of  commonly  made  approximations, 

•"  however,  the  radiative  transfer  calcula¬ 
tion  is  a  very  large  one.  In  fact  the 
radiative  transfer  calculation  is  more 
time  consuming  than  the  fluid  dynamic  cal¬ 
culation  by  a  substantial  factor.  This 
is  principally  due  to  the  larger  number  of 
independent  variables  required  to  describe 
the  radiation  field  as  compared  with  the 
fluid  field.  In  addition  to  spatial  co¬ 
ordinates  there  are  angular  and  frequency 
variables  as  well.  Consequently,  ill  cal¬ 
culations  of  the  time  -  dependent  behavior 
of  an  atmosnheric  fluid  dynamic  system 
have  employed  additional  simplifying  ap¬ 
proximations  in  order  to  reduce  the  amount 
of  calculation  to  solve  the  radiative 
transfer  equations.  In  our  presentation 
we  plan  to  report  the  tests  of  one  approx¬ 
imation  scheme  by  comparing  the  results  of 
our  benchmark  radiative  transfer  calcula¬ 
tion  with  that  of  the  Mintz-Arakawa  global 
circulation  model  as  used  at  the  RAND 
Corporation  for  climate  dynamics  calcula¬ 
tions. 

In  the  perspective  of  the  above 
discussion  the  objectives  of  our  investi¬ 
gation  are  as  follows: 

(1)  Formulate  a  theoretical  treat¬ 
ment  of  atmospheric  radiative  transfer 
oriented  to  the  determination  of  heating 
rates  in  the  atmosphere.  The  treatment 
must  be  applicable  to  arbitrary  conditions 
of  insolation,  surface  boundary  condition, 
and  atmospheric  distribution  of  gaseous 
and  particulate  absorbers  and  scatterers; 

(2)  Within  basic  well- justified 
approximations,  develop  a  computer  code 
which  incorporates  the  best  available 
techniques  and  data.  The  resulting  com¬ 
puter  code  should  qualify  as  a  standard 
for  comparison  with  radiation  routines 
contained  in  lluid  dynamics  codes.  Econ¬ 
omy  of  calculation  is  considered  to  be  of 
secondary  importance  compared  to  the 

JHabo  ve  objective; 
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(3)  Apply  the  computer  code  to  the 
calibration  of  a  radiation  routine  of  a 
climate  dynamics  general  circulation 
model  (GCM) . 

As  suggested  above,  there  are  sev¬ 
eral  different  classes  of  radiation  treat¬ 
ments.  The  radiation  subroutines  of  the 
fluid  dynamic  codes  contain  simplifying 
assumptions  in  order  to  achieve  computa¬ 
tional  speed.  These  calculations  must  be 
performed  repeatedly  for  many  spatial 
positions  and  at  many  time  steps  during  a 
typical  fluid  dynamic  calculation.  Ex¬ 
amples  of  the  treatment  of  radiation  for 
GCM's  are  discussed  by  Manabe  and  Stick¬ 
ler  (1964),  Sasamori  (1968)  and  Joseph 
(1966).  Radiation  in  the  planetary 
boundary  layers  is  considered  by  Atwater 
(1971).  At  the  other  extreme,  atmospheric 
spectroscopic  investigations  call  for  a 
very  detailed  treatment  of  narrow  fre¬ 
quency  intervals.  In  this  case,  interest 
in  angular  dependence  and  state  of  polari¬ 
zation,  in  addition  to  absorption  line 
profile  information,  produces  more  infor¬ 
mation  than  is  needed  for  atmospheric 
heating  and  greatly  increases  the  compu¬ 
tational  expense.  Examples  cf  this  ap¬ 
proach  are  to  be  found  in  the  extensive 
calculations  of  Dave  (1968),  Sekera 
(1963),  and  of  Thompson  and  Wells  (1971). 

Only  a  few  calculations  have  been 
leported  which  give  a  detailed  treatment 
of  the  entire  spectrum  with  the  objective 
of  evaluating  the  radiative  heating  of  the 
atmosphere.  Such  a  treatment  is  that  of 
Rasool  and  Schneider  (1971)  who  examine 
the  effects  on  global  climate  of  increases 
m  aerosols  and  CO2.  They  employ  high 
resolution  in  the  specification  of  the 
model  atmosphere  and  in  the  frequency  de¬ 
pendence  in  the  IR.  However,  the  angular 
distribution  and  interaction  of  scatter¬ 
ing  and  absorption  are  treated  very  ap¬ 
proximately,  Another  detailed  radiation 
calculation  of  planetary  atmospheres  is 
that,  of  Hunt,  and  Grant  (1569)  who  inves- 
tigate  the  effect  of  high  clouds  on  IR* 
radiation  emerging  from  the  atmosphere. 

They  calculate  accurately  the  radiative 
transfer  due  to  scattering  and  absorption 
but  use  an  idealized  phase  function  and 
low  frequency  resolution.  As  is  dis¬ 
cussed  in  greater  detail  below,  our  for- 
mulation  attempts  to  remedy  the  weaknesses 
of  botn  of  these  formulations. 

2.  FORMULATION  OF  THE  RADIATIVE  TRANS¬ 
FER  PROBLEM 

The  objectives  of  the  radiation 
code  discussed  in  the  Introduction  dic¬ 
tate  that  certain  approximations  which 
reduce  its  complexity  can  be  made.  We 
identify  these  approximations  in  this  sec¬ 
tion  and  comment  on  each  in  a  sentence. 

To  reiterate  the  objectives  of  the  study 
we  desire  an  accurate  treatment  of  the  ’ 
radiative  heating  rate  of  the  troposphere 
consistent  with  the  level  of  atmospheric 
data  obtainable  from  GCM  calculations  and 
global  observations.  These  data  are 


limited  at  present  but  we  anticipate  sub¬ 
stantial  improvement  in  them  in  the 
future  and  design  our  calculation  accord¬ 
ingly. 

Several  approximations  of  physical 
processes  have  been  made: 

(1)  The  emissivity  and  absorption 
coefficient  of  the  atmospheric  constitu¬ 
ents  are  assumed  to  be  that  appropriate 
to  local  thermodynamic  equilibrium.  This 
approximation ,  which  requires  molecular 
collisions  to  maintain  thermodynamic 
equilibrium  populations  of  those  states 
Pa t 1 i c ipa ting  in  radiative  transitions, 
is  valid  in  the  troposphere  where  the 
atmospheric  pressure  is  high. 

-(2)  The  adiation  field  can  be 
considered  to  be  unpolarized.  While 
molecular  and  Mie  single  scattering  events 
result  in  partially  polarized  radiation, 
it  has  been  shown  (Howell,  1970)  that  the 
error  resulting  from  a  polarization- 
independent  treatment  is  in  many  cases 
less  than  one  percent. 

(3)  The  air  is  assumed  to  be  non- 
refractive.  Refraction  effects  affect 
heating  rates  only  for  near-grazing 
angles  of  incidence  of  the  solar  beam  and 
these  angles  occupy  a  very  small  fraction 
climatological  events. 

(4)  Aerosols  are  assumed  to  be 
composed  of  uniform  spherical  particles 
with  a  homogeneous  index  of  refraction. 

It  is  known  that  this  assumption  is  not 
entirely  valid  for  atmospheric  aerosols; 
limited  investigation1  indicate  that 
scattering  is  not  stsongly  dependent  on 
shape  and  the  theory  of  non- spherical 
particles  is  very  much  more  difficult 
than  for  spheres. 

Geometrical  approximations  which 
permit  the  radiation  calculation  to  be 
treated  as  a  stratified  plane-parallel 
atmosphere  are  also  implied: 

(1)  The  atmosphere  and  boundary 
are  taken  to  be  plane-parallel.  The  at¬ 
mosphere  itself  is  well-approximated  as 
plane  when  tropospheric  radiation  (where 
mean  free  paths  are  short)  is  being  con¬ 
sidered,  when  the  clouds  are  stratified, 
and  the  sun  angle  is  not  near  grazing 
incidence  where  curvature  effects  should 
be  considered.  The  lower  boundary  must 
also  be  a  horizontal  plane. 

(2)  In  order  to  carry  out  an 
average  over  the  angle  of  azimuth,  it  is 
necessary  that  the  quantity  of  interest 
be  a  scalar ,  such  as  t lie  atmospheric 
heating  rate. 

A  number  of  mathematical  approxi¬ 
mations  have  also  been  made.  The  dis¬ 
cussion  of  these  (discretization  of 
angle,  frequency,  space)  takes  place  in 
the  next  section  along  with  the  descrip¬ 
tion  of  the  computer  code  itself. 
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3. 


uiisciurnoN  of  the  code 

The  code  has  several  unique  fea¬ 
tures,  amonf,  which  arc  flexibility  with 
respect  to  vertical  zoning,  ability  to 
attain  arbitrary  accuracy  by  adjusting 
a  few  parameters,  treatment  of  the  com¬ 
bined  line  absorption-scattering  problem 
and  particularly  of  band  overlap  regions, 
and  the  generality  of  its  surface  boundary 
condition.  These  features,  and  others, 
arc  described  in  more  detail  below. 

3 . 1  Transport  Formulation 


The  central  part  of  the  code  is 
that  which  solves  the  radiative  transfer 
equation.  The  algorithm  of  Grant  and 
Hunt  (1969)  has  been  chosen  for  this 
purpose.  This  algorithm  is,  in  our  esti¬ 
mation,  the  best  currently  available  for 
tli e  solution  of  the  general  monochromatic 
absorbing-scattering  problem.  It  el 
nates  altogether  the  need  for  a  cumb  i 
some,  and,  in  the  presence  of  clouds, 
prohibitively  expensive,  scattering 
iteration  procedure.  It  is  computation¬ 
ally  stable  and  economical,  allows  zones 
of  arbitrary  size,  conserves  flux  pro¬ 
vided  only  that  the  phase  function  is 
properly  "normalized",  and  guarantees 
positive  intensities  regardless  of  the 
optical  depth.  It  provides  for  error 
reduction  by  the  adjustment  of  a  few 
simple  parameters,  and  furnishes  error 
estimates  for  various  quantities  in¬ 
volved  in  the  computation. 

At  present,  the  code  assumes  that 
zones  are  homogeneous  with  respect  to  the 
albedo  for  single  scattoring  and  the 
phase  function;  this  permits  the  doubling 
method  to  be  used  for  reflection  and 
transmission  matrices.  However,  since 
neither  the  Planck  source  term  nor  the 
solar  source  term  (obtained  after  split¬ 
ting  the  intensity  into  a  solar  and  a 
diffuse  part)  are  constant  across  zones, 
the  normal  doubling  equations  for  source 
vectors  cannot  be  used.  We  have  derived 
doubling  formulae  for  these  inhomogeneous 
source  vectors  which  permit  them  to  be 
calculated  with  a  minimum  of  computing 
time . 


In' the  following  di scussi on  ,  space 
limitations  force  us  to  assume  that  the 
reader  is  familiar  with  the  Grant  and 
Hunt  (1969)  technique.  Consider  the  zone 
of  Figure  1,  composed  of  2*^  primary 
layers  of  optical  thickness  At.  Presume 
that  the  source  in  question  varies  across 
the  zone  according  to  f(t) ,  and  that  the 
source  vectors  for  pi imary  layer  i  are 


f(q) 


(i) 


where  the  Ep  are  independent  of  t,  and  £ . 
is  the  mid-point  cf  the  primary  layer,  1 


£.  *=  to  +  i At  -  ^6t  .  (2) 


Zone 


Fig.  1.  Zone,  homogeneous  with  respect 
to  single-scattering  albedo  and  phase 
function,  and  composed  of  2N  primary 
layers . 


Because  the  layer  superposition  formulas 
for  source  vectors  are  linear  in  the 
source  vectors,  the  superposition  of  any 
adjacent  2n  primary  layers  within  the 
zone,  beginning  at  the  ith  one,  yields 

2n 

£(i,n)  ’  E  Vi!,*  £<Wl>  « 

1=1 

(the  expressions  for  the  V's  are+not  of 
interest  here).  The  notation  ET.  .  is 
a  shorthand  for  ^ 1 » n J 

Ei  +  (iH)  +  ...  +  ( i  +  2n - 1 )  * 

For  the  solar  source,  we  have 

-t/u 

f(x)  =  e 


where  p  is  the  cosine  of  the  solar 
zenith  Sngle.  Assume  that  the  first  2n 
primary+laycrs  have  been  composed,  ro 
that  E7.  .  are  known.  The  next  2n 

primary''  *  J  layers  will  then,  by  Eqs. 
(2)  and  (3) ,  have  source  vectors 


+  *  +  -[r+U+2  ^t-^AtJ/p 

,n  ,  =  >  V"  e  0  1  0 

(1+2  ,n)  /  v  n,l 

1=1 


h  E  7  x  x 
n  (1  ,n) 


(4) 


whe  re 


hn  -  e 


- 2nAr/ p 


(5) 


Using  F.q .  (4),  the  source  vectors  for  the 
combined  layer  of  2n'*1  primary  layers  arc 
given  by  the  (source  superposition)  for¬ 
mulas  <“ 
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E(l,n+1)  =  tnrn[rn  E(l+2n,n)  +  E(l,n)] 


+  S(i*2n,n) 


more  tractable  function  than  Bv[T(t))  is 
needed.  For  zones  across  which  t lie  tem¬ 
perature  change  AT  is  nr,?  large  (AT<l<i°K) 
it  can  be  shown  by  expanding  the  Planck 
function  that 


Vn[hnrn  E(l,n)  +  E(l,n)] 


+  hn  E(l,n; 


r(l,n+l)  ‘  tnFn  [r n  E(l,n)  +  E(l+2n,n)] 


(l,n) 


Vn[rn  Ea,n)  +  hn  E"(l .»)] 


'(1  »n) 


where  rn  and  tn  are  the  reflection  and 
transmission  matrices  for  any  layer  of 
thickness  2n  At  within  the  zone,  and 


r  =  (I  -  r  r  ) 
n  v  n  nJ 


From  the  definition  (5)  of  h  , 

n 

Vl-hn  • 


Together,  F.qs,  (6-8)  form  a  doubling 
scheme  for  the  solar  (or-  any  exponential- 
in-r)  source.  The  scheme  is  initialized 
by 


1(t)  =  B  +  B'(t-t  ) 

0  .  .  0 

reasonably  approximates  the  variation  of 
across  the  zone  of  Figure  1,  where 

\  B  VT<V]  «  and 


B  [T(t  ))  -  B  [T(t  )] 
r  i  =  _ i _ v  o 

0  T  -  T 

1  0 

Consider  a  sublayer  of  thickness  2n  At 
(0<n<N)'  beginning  at  t=t0.  Across  this 
sublayer,  write  f(r)  as  a  linear  combi¬ 
nation  of  a  constant  part  and  an  anti¬ 
symmetric  part: 

f(t)  =  anf"  +  B'f"  ,  t  <  t  <,  t  +  2n  At 


where  f 


a  ■=  B  +  2n_1  At  B' 
n  o 


and  f£  =  t  -  2n  E  At 


We  now  derive  the  doubling  formulae  for 
the  source  vectors  1'^  and  Z,*  correspond¬ 
ing  to  f"  and  f“.  The  major  advantage  of 
working  with  Y  and  Z  is  that  they  are 
symmetric  and  anti -symmetric ,  respective¬ 
ly: 


Y  =  Y  =  Y 
n  n  n 


-  Z„  =  Z 
n  n 


— -F  (ill  ,u  ) 

V,  vv 


(1,0)  4n 


Svw  *(T  +7*T)/p 
7T-  e  0  1  0  At 


-At/u 
h  =  e  0 


• — P  (±p  ,p  ) 
Pm  v  m  « 


(these  relations  may  be  Moved  inductive¬ 
ly)  .  The  doubling  formula  for  Yn  is 
Straightforward,  since  the  corresponding 
source  fa  is  homogeneous  across  the  zone: 

Yn  +  1  "  (r  +1)  «  1 ) Y  .  (10) 

n+1  1  n  n  n  '  Jn  ' 

Consider  the  source  vector'  Z 
corresponding  to  the  source  f’,nctionn+i 

f?+E  =  T-2n  At-t  ,  t  <  t  <  t  +2n+E  At 

Li  a  n  *—  —  ft 


and  iterated  from  n-0  to  n=N-l.  (S  is 
the  solar  flux  at  frequency  v,_w  isvthe 
albedo  for  single  scattering,  Pv  is  the 
azimuthally-averaged  phase  function,  and 

P., . P  are  the  angular  quadrature 

points . ) 

For  the  Planck  source, 
f(T)  -  Bv[T(t)] 

where  Bv  is  the  Planck  function  and  T  is 
the  temperature.  In  order  to  derive 
doubling  formulas  for  this  source,  a 


We  treat  this  source  function  separately 
in  the  two  halves  of  the  interval: 

f£+E  =  f£  -  2n  E  At  ,  T( <t<t ^ +2n  At  , 

f?  +  E  =  T-2n  At-t  ,  t  +2n  At<t<t  +  2n*^  At 
l  o  o  —  —  o 


15  t'-t  b  f?  +  2n  E  At  , 
o  L 


IfK- 


T  <t'<t  +2  At  . 
0-  —  0 


L 


I 


The  correspond ing  source  vectors  are 
E0,n)  '  K  -  2"'1  4’  \  • 

E(l+Zn ,n)  '  ^  *  2"-1  *„  • 

Employing  the  source  vector  superposition 
formulae , 

Z_.i  =  t  T  [r  £,,.,n  s  +  £*  .1 

n+1  n  n 1  n  (1+2  ,n)  (l,n)J 


+  £d  +  2n,n) 


wherA 


»  [I+t„r  (I-r  )]Z„ 
1  n  n v  n' 1  n 


+  Sn^Vn^V^n 


On‘l  . 

2  At 


(ID 


Clearly, 


en+l 


=  2gr 


(12) 


Together,  Eqs .  (10-12)  provide  a  doubling 
scheme  for  the  Planck  (or  any  linear- in- t) 
source.  The  scheme  is  initialized  by 

1/P. 


V  =  (I-w)At 


Z  =  0 
o 

0  ** 


1/p. 


and  iterated  from  n-0  to  n=N-l,  The  final 
values,  Yjyj ,  Z^,  and  g^,  are  combined  just 
as  the  source  functions  in  Eq .  (9)  are 
combined: 


Z(1,N)  =  aN  YN  +  B’  ZN 


+  8m  B’)Yw  +  B'  Z, 


Z(1,N).  =  aN  YN  +  B'  ZN 


c  <B,  +  8N  B')Y„  -  b-  zn  . 

More  accurate  pieccwise-poly- 
nomial- in- t  approximations  to  Hv[T(t)J 
are  possible,  at  tlie  expense  of  more  com¬ 


plicated  source  doubling  formulas.  In 
particular,  we  have  derived  economical 
quadratic  and  cubic  schemes  for  future  in¬ 
clusion  in  the  code  (if  the  increased  ac¬ 
curacy  so  obtained  warrants  their  use) . 

The  solar  and  Planck  source  vec¬ 
tors  are  simply  added  to  give  the  total 
source  vector. 

3 . 2  Material  Properties 

The  properties  of  the  various 
gaseous,  liquid,  and  solid  constituents 
of  the  atmosphere  enter  the  calculation 
of  the  absorption  and  scattering  coeffi¬ 
cients  and  phase  function. 

3.2.1  _  Scattering  Treatments  —  Rayleigh 

scattering,  with  the  non-negligible  de¬ 
polarization  correction,  is  included  ol- 
lowing  Penndorf  (1957).  A  full  Mie 
scattering  calculation  (for  homogeneous 
spheres),  with  integration  over  a  user- 
specified  aerosol  size  distribution,  is 
also  included.  Any  of  the  commonly-used 
r  analytic  aerosol  size  distributions 

(modified  gamma,  Gaussian,  log-normal)  are 
available  as  options,  or  the  user  may 
supply  these  data  in  card  form.  The  com¬ 
plex  index  of  refraction  of  water  is  tabu¬ 
lated  over  the  entire  frequency  range  of 
interest  in  the  atmosphere,  so  that  the 
accuracy  with  which  the  most  important 
aerosols,  clouds,  are  treated  depends  en¬ 
tirely  on  the  accuracy  with  which  the 
size  distribution  is  known.  The  indices 
of  refraction  of  other  aerosol  substances 
are  tabulated  as  completely  as  available 
data  permit. 

Since  the  Mie  scattering  calcula¬ 
tion  can  be  costly,  especially  when  clouds 
are  present  and  the  phase  function  must  be 
calculated  at  many  angles,  an  option  is 
provided  to  use  the  Heny  y-Greenstein 
phase  function  instead.  The  Mie  phase 
function  is  then  calculated  at  a  rela¬ 
tively  small  number  of  angles  and  the 
single  parameter  in  the  Henyey- Greenstein 
function  chosen  to  give  a  best  fit  at 
those  angles.  Several-parameter  analytic 
phase  functions  are  being  sought  in  order 
to  give  greater  freedom  in  this  fitting 
process. 


The  Mie  scattering  sections  of  the 
code  were  checked  against  the  tables  of 
Deirmendjian  (1969).  In  the  integration 
over  particle  size,  Deirmendjian  picked 
integration  intervals  and  tabular  angles 
"by  eye."  In  order  to  perforin  this  inte¬ 
gration  by  computer  program  an  alterna¬ 
tive  procedure  was  devised.  It  involves  a 
crude  preliminary  size  integration,  which 
is  used  to  limit  the  size  of  the  total 
interval  of  integration,  to  pick  the  n- 
tegration  mesh,  and  to  establish  the  an¬ 
gular  mesh  on  which  the  phase  function  is 
to  be  calculated.  A  Romberg  integration 
is  then  done  to  obtain  the  final  results. 


If  the  phase  function  lias  too 
large  a  forward  peak,  it  is  truncated 


I 


after  the  fashion  of  Potter  (1970)  who 
fimls  that  the  resulting  intensity  er¬ 
rors  arc  generally  less  than  one  percent 
Mux  errors  arc  even  smaller.  By  thus 
ensuring  that  the  phase  function  does  not 
have  a  strong  angular  dependence,  accura¬ 
cy  can  he  obtained  in  fluxes  using  a  re¬ 
latively  coarse  angular  specification  of 
the  intensity  (3-7  Gaussian  angles  are 
currently  being  used) . 

3.2.2  Absorption  Treatment  -  In  addition 
to  the  continuum  absorption  due  to  Mie 
scattering,  two  other  absorption  continua 
are  included  in  the  code:  the  ll20  contin¬ 
uum  in  the  8-13p  window  region;  and  the  N-> 
continuum  around  4p,  Data  for  both  of 
these  continua  arc  taken  from  McClatchey, 
et.al.  (1970).  McClatchey 's  group  has 
also  furnished  what  we  find  to  be  the  best 
currently  available  transmission  functions 
for  molecular  bands  (barring  a  prohibi¬ 
tively  expensive  line-by-line  calcula¬ 
tion).  These  transmission  functions  are 
tabular  fits  to  actual  line-by-line  cal¬ 
culations,  with  avei aging  intervals  of 
20  cm  ,  (Belov  350  cm"^  we  use  the 
Goody  random  band  model  for  the  H20  rot; - 
tion  band.-)  Transmissions  arc  calculated 
separately  for  II20,  for  0,,  and  for  the 
uniformly  mixed  absorbers  (C02,  N20,  CO, 
CH4,  and  02).  The  scaling  approximation 
is  used  with  parameters  chosen  to  give 
the  best  agreement  with  line-by-line  cal¬ 
culations.  The  results  are  claimed  to  be 
accurate  to  better  than  10  percent,  and 
in  most  cases  to  better  than  1-2  percent; 
the  worst  errors  are  in  the  wings  of 
strong  bands. 

3 • 3  Splitting  Into  Monochromatic 

Problems 


M 


wu> a  E  *i 

i=l 


-k.u/p 

e 


(14) 


where  M  is  small,  0(10)  say  (thus  preclud¬ 
ing  the  sum  from  being  an  ordinary  quadra¬ 
ture  for  the  integral ,  which  would  require 
hundreds  or  thousands  of  terms) .  If  iv  is 
the  diffuse  intensity,  suppose  that  cor¬ 
responding  to  Eq.  (14)  we  have 

M 

aav  =  h  j£v  \  dv  s  Zai  ii  •  (1S) 

i=l 


Each  i^  is  chosen  to  satisfy  a  monochro¬ 
matic  problem  with  absorption  coefficient 


(where  2  is  vertical  distance).  Recombin¬ 
ing  the  ii  according  to  Eq .  (15),  with  the 
_  a*  determined  by  some  fitting  scheme  lead- 
•’  ing  to  Eq .  (14) ;  produces  the  frequency- 
averaged  intensity  ifi  . 


Yamamoto  (1971)  explicity  assumes 
and  Grant  and  Hunt  (1969)  tactily  assume 
fnat  the  above  procedure  is  correct.  How¬ 
ever,  they  do  not  actually  give  a  deriva¬ 
tion  beginning  from  the  radiative  transfer 
equation.  The  present  authors  find  that 
a  further  assumption  is  needed,  namely  that 


-kvu/p 

e 


dv 


M 


i*=l 


-kiu/p 


(16) 


The  Grant  and  Hunt  algorithm  is 
applicable  to  monochromatic  problems.  Un¬ 
fortunately,  in  solving  atmospheric  radia¬ 
tive  transfer  problems,  frequency  groups 
may  not  be  taken  narrow  enough  that  the 
absorption  coefficient  remains  reasonably 
constant  across  them.  Typically,  a  fre¬ 
quency  group  will  be  upwards  of  20  cm-1 
in  width  and  will  contain  numerous  ab¬ 
sorption  lines.  Therefore,  a  further  ap¬ 
proximation  is  necessary. 

The  approximation  we  have  chosen 
is  suggested  by  some  work  of  Grant  and 
Hunt  (1969)  and  of  Yamamoto  (1971).  The 
idea  expressed  in  Eq.  (14)  antedates  their 
work,  however,  and  may  be  traced  to,  for 
example,  Cowling  (1950). 

Let  the  transmission  function  be 


TAv<u>  "  Sv 


L 


/kvu/p 

e  dv 


(13) 


Granting  this,  the  splitting  into  mono¬ 
chromatic  problems  follows.  We  argue  for 
the  plausibility  of  Eq.  (16)  based  on  cer¬ 
tain  ideas  about  Lebesgue  quadrature  and 
about  the  behavior  of  iv .  This  point  is 
considered  more  definitively  by  Freeman, 
Wiscombe,  and  England  (1972). 

It  is  difficult  to  obtain  the  fit 
in  Eq.  (14)  in  an  automatic,  stable,  con¬ 
vergent  fashion,  since  the  problem  of  ex¬ 
ponential  -  sum-approximation  is  notorious 
for  its  numerical  ill-conditioning  and 
general  intractability  (cf.  Lanczos,  1956). 
Recently ,  however,  Cantor  and  Evans  (1970) 
devised  a  way  of  "factoring  out"  the  ill- 
conditioned  part  of  the  problem,  and  their 
method  has  been  incorporated  into  the  code. 
It  is  accurate,  computationally  fast,  has 
guaranteed  convergence,  and  provides  the 
unique  best  fit  in  the  least  -  squares  sense. 

3.4  Boundary  Conditions 


,  .  The  surface  boundary  condition  is 

wheie  u  is  the  scaled  absorber  amount.  -4  on  formulated  in  complete  generality,  allow- 

1-lt  TAv  Wlt!l  a  SUI"  °f  exponentials,  X*7o^ng  for  the  specification  of  the  full  bi¬ 

directional  reflectivity.  There  are  also 
provisions,  in  case  less  information  is 


available,  to  use  the  directional- 
hemi spherical  reflectivity  or  the 
hemispherical  reflectivity  (albedo).  In 
the  latter  two  cases,  diffuse  reflection 
is  assumed.  The  emissivity  is  calculated 
from  the  reflectivity  using  Kirchoff's 
Law.  In  case  only  the  directional  or 
hemispherical  emissivity  is  available,  the 
directional-hemispherical  reflectivity  or 
albedo  is  calculated  therefrom  using 
Kirchoff's  Law. 

Theoretical  derivations  of  rough- 
surface  bidirectional  reflectivities,  with 
roughness  specified  by  Cox  and  Munk's 
(1956)  distributions  for  the  sea  surface 
and  by  McStravick's  (1972)  model  for  land 
surfaces,  have  been  incorporated  into  the 
code.  Various  albedo  data,  such  as  that 
of  Knnov  (1947),  have  also  been  in¬ 
cluded  . 

4.  TEST  PROBLEMS 

There  arc  very  few  exact  solutions 
in  ra'diative  transfer  theory.  Those  that 
do  exist,  are  for  monochromatic  problems 
only.  We  have  tested  the  code  against 
several  such  problems.  For  the  case  of 
no  scattering,  and  using  a  realistic  ab¬ 
sorption  coefficient  (a  Lorcntz  line)  with 
pressure  and  temperature  dependence,  prob¬ 
lems  were  run  with  various  temperature 
profiles  and  surface  boundary  conditions 
and  compared  with  the  analytic  solution. 

For  a  scattering  test  problem,  comparisons 
were  made  with  the  Rayleigh  scattering 
solutions  tabulated  by  Sckcra  and  Kahlc 
(1966).  Agreement  in  these  cases  was  ex¬ 
cellent. 

Investigation  is  continuing  into 
the  accuracy  of  the  approximation  of  Eq. 
(16).  A  detailed  linc-by-linc  transfer 
calculation  is  being  designed,  using  a 
model  line  structure,  to  ascertain  the 
limits  of  its  validity. 

*5.  COMPARISON  WITH  GCM  CALCULATIONS 

An  application  of  our  code  is  being 
made  to  the  calibration  of  the  radiation 
subroutine  of  the  Mintz -Arakawa  GCM  as 
used  in  the  climate  dynamics  research  pro¬ 
gram  at  tlie  Rand  Corporation  (Gates, 

1971) .  Results  of  comparisons  between 
the  radiation  treatment  in  that  code  and 
our  code  will  be  included  in  our  presen¬ 
tation.  These  will  consist  of  compari¬ 
sons  of  radiative  heating  rates  of  the 
atmosphere  for  selected  locations  on  the 
global  surface.  Atmospheric  input  data 
arc  taken  from  the  GCM  calculation. 

6,  REFERENCES 

Atwater,  M. ,  1971:  The  radiation  budget 

for  polluted  layers  of  the  urban 

environment.  J.  Appl.  Meteor..  JO 

205-214.  ’  ~ ’ 


Cantor,  1).  and  J,  Evans,  1970:  On  approx 
imation  by  positive  sums  of  powers 
SIAM  J.  Appl.  Math.,  L8,  380. 

Cowl ing ,  T, G . ,  1950:  Atmospheric  absorp¬ 
tion  of  heat  radiation  by  water 
vapor.  Phil.  Mag.,  £1,  109. 

Cox,  C.  and  W.  Munk,  1956:  Slopes  of  the 
sea  surface  deduced  from  photo¬ 
graphs  of  sun  glitter.  Bull. 
Scripps  Inst.  Oceanog.,  6,  401. 

Have,  J.V.,  1968:  Transfer  of  solar 

ultraviolet  radiation  through  the 
earth's  molecular  atmosphere. 
JQSRT,  8,  25-38. 

Deirmendjian,  D.  ,  1969:,  Electromagnetic 
scattering  on  spherical  polydis- 
persions.  Elsevier. 

Freeman,  B.E.,  W.J.  Wiscombe,  and  W.G. 
England,  1972:  The  effects  of 
meso-scale  and  small-scale  inter¬ 
actions  on  globil  climate.  Sys¬ 
tems,  Science  and  Software  Report 
3SR-1034,  La  Jolla,  Calif. 

Gates,  W.L.,  E.S.  Batten,  A. B .  Kahle,  and 
A.B.  Nelson,  1971:  A  documenta¬ 
tion  of  the  Mintz-Arakawa  two- 
level  atmospheric  general  circula¬ 
tion  model.  Advanced  Research 
Projects  Agency  Report  R-877. 

Grant,  1.  and  G.  Hunt,  1969:  Discrete 

space  theory  of  radiative  transfer 

I.  Fundamentals.  Proc.  Roy.  Soc. 
Lond.  A,  3£3,  183. 

Howell,  H.  and  H.  Jacobowitz,  1970: 

Matrix  method  applied  to  multiple 
scattering  of  polarized  light. 

J.  Atmos.  Sci.,  £7,  1195-1206. 

Hunt,  G.  and  I.  Grant,  1969:  Discrete 

space  theory  of  radiative , transfer 
and  its  application  to  problems  of 
planetary  atmospheres .  J.  Atmos. 
Sci.,  26_,  963-972. 

Joseph,  J.,  1966:  Calculation  of'  radia¬ 
tive  heating  in  numerical  general 
circulation  models.  Dept,  of 
Meteorology  Tech.  Report  No.  1, 
University  of  California,  Los 
Angeles . 

Krinov,  E.L.,  1947:  Spectral  reflectance 
properties  of  natural  formations. 
Nat.  Res.  Council  of  Canada  Tech. 
Trans.  439,  Ottawa. 

Lanczos,c.f  1956:  Applied  Analysis. 

Prenti  ce-llall ,  p.  272. 

Manabc ,  S.  and  R.F.  Strickler,  1964: 

Thermal  equilibrium  of  the  atmos- 
1  Q/J  phere  with  a  convective  adjustment. 
J.  Atmos.  Sci.,  21  ,  361-385. 


McClatchcy,  R.,  ct.a].,  1970:  Optical 
properties  of  the  atmosphere. 

Air  Force  Cambridge  Report 
AFCKL-70-0527 ,  Cambridge,  Mass. 

McStravick ,  D.M, ,  1972;  The  effect  of 

surface  roughness  on  the  reflected 
and  emitted  energy  from  a  rough 
surface.  Pli.D.  Thesis,  Rice  Univ. 
Texas . 

Pcnndorfj  R,,  19S7:  Tables  of  the  refrac¬ 
tive  index  for  standard  air  and 
the  Rayleigh  scattering  coefficient 
for  the  spectral  region  between  0.2 
and  20p.  J.  Opt.  Soc.  Amer. ,  47, 
176.  — - 

Potter,  J.,  1970:  The  delta  function  ap¬ 
proximation  in  radiative  transfer 
theory,  J.  Atmos.  Sci.,  17_,  943. 

Rasool,  S.I.  and  S.H.  Schneider,  1971: 
Atmospheric  carbon  dioxide  and 
aerosols:  Effects  of  large  in¬ 
creases  on  global  climate.  Science, 
173,  138-141. 


Sasamori,  T.,  1968;  The  radiative  cool¬ 
ing  calculation  for  application  to 
general  circulation  experiments. 

J.  Appl .  Meteor.,  721-  729. 

Sckcra,  Z.,  1963.  Radiative  transfer  in 
a  planetary  atmosphere  with  imper¬ 
fect  scattering.  The  Rand  Corp. 
Report  R-1413-RR,  Santa  Monica, 
Calif. 

Sckera,  Z.  and  A.  Kahlc,  1966:  Scatter¬ 
ing  functions  for  Rayleigh  atmos¬ 
pheres  of  arbitrary  thickness. 

The  Rand  Corp.  Report  R-452-PR, 
Santa  Monica,  Calif. 

Thompson,  B.C.  and  M.B.  Wells,  1971: 

Scattered  and  reflected  light  in¬ 
tensities  above  the  atmosphere. 
Appl.  Optics,  1£,  1539-1549. 

Yamamoto,  G.,  M.  Tanaka,  and  S.  Asano, 

1971:  Radiative  heat  transfer  in 
water  clouds  by  infrared  radiation. 
JQSRT,  11_,  697. 


195 


